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The study of the static configuration of capillary surfaces as 
well as their dynamic behaviour has wide applications in many 
technological areas especially in space technology. The capillary 
surface is governed by the following "Laplace-Young" equation 
[Myskis A.D. : Low Gravity Fluid Mechanics (Spr inger-Ver lag) 

1987.3 


4^1 • 

I ■^l + |Vu|^ J 


ixyu = c i n O 


<P-1) 


Sfu.n I 

Z— = cos 0 on ^ 

■^l + |Vup 

When Q is known a-priori the problem amounts to the solution of a 
system of non-linear equations. In view of the early development 
there has been an extensive study even in standard monograph form 
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dealing with the question of the existence of solution of (P-l> 
[Finn R. : Equilibrium Capillary Surfaces (Spr inger-Ver lag) 1986]. 

The present thesis concerns itself with the computational 
aspects of the solution of <P-1> and eigen-analysis of 
oscillations of an ideal liquid contained in rigid containers 
under low-gravity conditions. 

(1) Generally O is not known a priori and therefore the problem 
becomes a "free-boundary" problem. In such a case the problem of 
the computation of the capillary surface with volume constraint on 
the liquid becomes quite complicated. In the thesis the above 
problem is treated invoking shape-optimization formalism and a 
mathematical justification of the method is derived. This seems 
to be the first attempt to solve the problem in this manner. 

The basic idea behind the application of this formalism 
consists of the following: 

For an admissible set H determine Ci € M s.t. 

J<0*> < Jidy V Q € M 

Where -7(0) is the energy functional for O. Here Q corresponds to 
a shape of the domain and it is varied over an admissible set M. 

We derive a mathematically justified numerical procedure based on 
ideas outlined in the work of Pironneau CPironneau 0 : Optimal 

Shape Design for Elliptic Systems. (Spr inger-Verlag) 1984]. 

(2) Earlier numerical simulation studies using finite-elements 
have been carried out for <P-1) when O is known a priori [Brown 
R.A. .-[Finite - Element Methods for the calculation of Capillary 



V 


Surfaces. Jour of Comp. Phys. 35 (1979)] 

In this situation also, inspite of the presence of high gradients 
in certain regions no adaptive scheme was available. Hence 
possibly for the first time a reliable adaptive strategy for the 
/i-version finite-elements has been formulated. It is based on 
some error indicators suitable for this problem in particular. 

(5) Then we consider the problem related to the study of 
oscillations of an ideal liquid in a rigid container. The 
equation governing the liquid oscillations are as follows (Myskis 
A.D. et al.: Low-Gravity Fluid Mechanics. (Springer- Verlag) 1987]. 
=0 in O 

If ° ® E 

’^( If] * (^] * ^ *" '' 

l?[lf]^- If =" -- 

= 0 [volume constraint condition] 

There has been quite an amount of work in this direction using 
semi -analyt i c methods [Bauer H.F. : Linear Liquid Oscillations in 
Cylindrical Container under Zero-Gravity. Technical Report. 
Universit&t Bundeswehr MUnchen (1989)]. 

However the range of applicability of these methods is very 
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limited. Therefore we propose a finite-element method for the 
calculation of eigen-values based on the variational formulation. 
The most distinguishing feature of the above problem is the 
presence of higher-order differential operator on the boundary 
surface. Therefore it does not allow standard procedure to reduce 
it to a variational form. However this has been achieved by 
introducing two distinct operators which together help in reducing 
the problem to a standard eigen-value problem. 

<4) The problem of computation of eigen-values for vibration of 
liquid in a container containing a large number of rigid tubes is 
very important form the point of applications. We deal with this 
problem using asymptotic methods based on homogenization 
principle. The problem is fairly complicated due to its 
5-dimensional nature. For solving this problem we make use of 
asymptotic analysis in the spirit of Lions CLions J.L ; Some 
Methods in the Mathematical Analysis of Systems and their Control 
(Gordon Breach) 1981] and construct the homogenized operators. 
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CHAPTER I 
INTRODUCTION 

1.1 MOTIVATION ! 

With the development of space-technology and the 
diversification of investigations conducted , the role of the 
fluid is becoming very important. The space-vehicles contain 
liquid fuel tanks and many other liquid containing devices. 
Various types of fluids are used in the rocket engines power 
units, life support systems, in the temperature control 

equipments etc. The static fluid configurations can provide 
essential information for storage. capillary pumping, venting 
etc. whereas the dynamic analysis of a liquid is very important 
for the design of the structural systems and more importantly for 
the stability of the system. Hence, the design and operation of 
such systems calls for the theoretical and experimental 
investigations of the behaviour of the liquid under low-gravity 
conditions. However, conducting experiments in space is expensive 
whereas the simulation of zero-gravity under terrestrial 
conditions is complicated and expensive. Thus there is a growing 
need for a rigorous mathematical analysis along with computational 
st-udies in order to produce substantial results so that the fluid 
behaviour may be predicted under extreme conditions. 

There is a special interest to study the static and 
dynamic behaviour of the liquids in a vessel under low gravity 
conditions. The present study has been aimed at liquids in rigid 
containers and shall not deal with problems of liquid drops, 
bubbles, liquid bridges etc. which may be important from the 
practical point of view. As the natural frequencies of 
oscillations of the liquid are much lower than the structural 
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frequenci.8 the vessel is considered to be rigid. 

The basic problem concerning the static state is the 
determination of the capillary surface of equilibrium. The 
capillary surface is governed by the following "Laplace-Young" 
equat i on . 


V XL 


l + jVup 
V .n 


+ i>o c c in O 


(P-1 > 


l+lVu| 


cost © on aCi 

C 


where m stands for the surface, ^>o is a system parameter 
depending on the gravity. the angle of contact (for details 

refer to Chapter IV> 

If H denotes the mean curvature of the surface 

CNitsche (1976). Myskis (1987)3, then 2« = V f " "" 1 

L ■^l + IVul^ -I 


Then it can be seen that IH = c - hoM 
is obtained from (P-i>. 


( 1 . 1 ) 

Since the mean curvature can be 
explicitly expressed in terms of a function of the space 
variables, the capillary surface is often described as a surface 
of prescribed mean curvature with a given contact angle on the 
boundary. In the case when the gravity is zero one has bo = O and 
one can find from (1.1> that N is constant. Thus the capillary 
surface is termed as a surface of constant mean curvature or the 
so-called M surfaces. If addition there is no volume 
constraint Cie. C = O 3 one arrives at the problem of minimal 
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surfaces or the Plateau's problem. The problem in its generality 
is to solve (P-l> when Ci is not known a priori. This is termed in 
literature as the "free-boundary" problem whose numerical solution 
is in general difficult to obtain. Possibly the first published 
attempt to numerically solve the 3-D meniscus is by Orr et al 
[19773. They adopted a direct method to solve <P-1) using 
finite-element method. The above work was essentially a 
simulation study without any attempt to justify the method. Under 
these circumstances there is a need to develop a numerical scheme 
with a proper mathematical basis and to show the practical 
viability of the method. 

In earlier numerical simulation studies finite-element 
method has been applied by some authors to solve<P-l> when Q is 
known a priori. In view of the presence of high gradients of v 
in certain regions of liquid container domain there is a need for 
the improvement in the methodology. In other words, one has to 
change the model, mesh or the approximating structure of the 
computational methods so as to improve the quality of the 
solution. This is usually termed as adaptive refinement. 
Unfortunately, there were no adaptive schemes available for the 
problem. Thus, there is a need to develop a practical adaptive 
scheme for the problem. In this connection it is to be noted that 
as the operator occuring in <P-1> is non-linear, it is a 
non-trival task to obtain a reliable criteria for adaptive 
refinement. The study of the oscillations of the liquid is 
complicated from theoretical and computational point of view. The 
equations governing the oscillations of an ideal liquid are as 


follows : 



4 



where ^ stands for the velocity potential of the fluid. 
CO the square root of the eigen value. C a constant corresponding 
to the Lagrangian multiplier for the volume constraint, n the 
outer normal and ® the tangent vector to the surface at the line 
of contact and X a system parameter defined on the line of 
contact. <For details refer to Chapter V> . 

Earlier some computational studies have been carried out by Bauer 
C1965a. 1965b. 1981. 1981. 1989a. 1989b. 1989c] using 
semi -anal yt i c methods. The method is not suitable for complex 
geometries or even problems of higher dimensions. For some Soviet 
works reference can be made to Myskis C1987]. The finite - element 
method based on variational formulation can serve as an effective 
method to overcome these difficulties. One basic difficulty in 
the application of this method arises from the fact that the 
differential operator appearing on the boundary r is of higher 
order than the differential operator in the interior of the domain 
(O) . Thus the present structure of (P-2> is not readily amenable 
to a simple application of finite element method. Thus there is a 
necessity to develop a procedure in order to reduce <P-2) to a 



standard eigen-value problem in the variational formulation. 

The computation of the eigen-frequencies of the liquid 
oscillations containing tube bundles has important applications. 
They are useful in the study of reactors where the number of tubes 
is of the order of several tens of thousands. A substantial 
amount of work has been 'devoted to the computational aspects of 
the study of vibration of the tube-bundles by Planchard [1980. 
1982, 198?]. Planchard [1980,198?] applied the theory of 

homogenization to compute the eigen-frequencies. In low-gravity 
environment the capillary surface oscillations become important as 
the surface-tension assumes a dominating role. The problem is 
?-dimens ional in nature and homogenization is not a simple task 
requiring the construction some special operators. 

1.2 NATURE AND OUTLINE OF THE WORK 

Keeping in view of the fact that the capillary problem 
is an old one and some amount of work has been done, it has been 
proposed to address the following four different aspects of the 
capillary surface problem. 

1) The problem of construction of the capillary surface using 
shape optimization technique . 

2) The problem of adaptive refinement scheme for the Finite 
Element Method . 

?) The problem of the computation of the eigen-frequencies of 
oscillations of an ideal liquid under low-gravity . 

4) The problem of the computation of the eigen-frequencies of an 
ideal liquid containing a rigid tube-bundle using the theory 


of homogenization. 
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The above problems are presented along with the state of 
the art. methodology, implementation and results in chapter III. 
IV. V and VI repectively. In addition. the mathematical 
preliminaries are given in chapter II and overall conclusions in 
chapter VI I . 

'The present work develops numerical schemes giving proper 
mathematical justification. Adequate numerical computations have 
been carried out to support the practical viability of the 
methods. The software has been developed in FORTRAN language and 
the computations have been carried out on H.P. 9000/850 Mainframe 
computer and Graphic workstations. 

The contents of Chapter III and IV have been accepted 
for publication in the Journal "Numerische Mathematik" and Journal 
"Computer Methods in Applied Mechanics and Engineering" 
respectively, under the following titles : 

1. "A Shape-Optimization Technique for the Capillary Surface 
Problem" . 

2. "A Finite-Element Adaptive Strategy for the Capillary Surface 

Problem. " 

The contents of Chapter V have been submitted to Journal 
"Computer Methods in Applied Mechanics and Engineering" 
under the following title. 

5. "A Numerical Method (using Finite-Elements) for the 
Computation of the Eigen-frequencies of Oscillations of an 
Ideal Liquid in Low Gravity." 



CHAPTER II 


MATHEMATICAL PRELIMINARIES 

The aim of this section is to present some important 
classical results and definitions which will be useful in the 
seque 1 . 

^ Tit 

We shall denote n as an open set in [R with boundary Xi. 
Then we have the following definitions : 

C <0) denotes the family of functions with domain Q whose m^ 
derivative is continuous. 

C®*<0) denotes the family of infinitely differentiable functions. 

D<0) denotes the sub-set of C®*(Q> which have compact support in O. 

We denote Z.^(0> for 1 < /> < eo as the space of equivalence classes 

of functions s.t. f\/\^ ^ • I The symbol J <.) shall be 

used often instead of the complete notation for 

simplification purposes and the integration is understood to be in 

the Lebesgue sense] 

00 

When p = 00 , L <Q) denotes the space of equivalence classes of 

functicns which are bounded almost everywhere <a.e.>. It is a 

Banach space with the norm defined by |/| = srup J/Cx^l 

X e Q 

For p, <? > 1 and — + ^ = 1 we have the Hblder inequality 


^ |/IlP I'li'’ 


V / « L'’<n> and d«t‘’<Q) . 


The multi-index a is defined as the n-tuple 


a = (a,. a-, o ) 0 ot. are integers 

JL c n 1 
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Associated with the multi-index we have the following symbols 


l«l 


a , 




a 


Xjl. 


a 

. .X n 
n 


K 




We say that two mu 1 1 i - i ndi ces a and ft are related hy a. < ft i f - 
V 1 < 1 < n. 

^ i 

Finally we express the differential operator through the 
mul t i - index . 




dhl 

ax, I . . . ax n 
i n 


With the help of the above notations we can now define the Sobolev 

space 1*^<0) as the following 
P 


= fu € L^CCi> 
P I 


c L^ca> V 



equipped with the norm 


hl».P.o = C 

|o(|^ 

Vh»n p=2 the space will be denoted by 

and for any n € W*^(0) we denote its norm by i.e. 

An extension operator for l/^<0) is a bounded linear operator 

n P 
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Op p 

s.t = w V ^e e l*^<Q) 

O ifi p 

The norm of will depend on O in general. Next we shall define 
a class T of domains for which the norms of corresponding 
extensions do not depend on O c T. In such a case we say that O € 
T has the so-called uniform extension property. 

This extension property is summed up in the following theorem due 
to Chenais C1975] 

Theorem 2.1 - 

Let 0,h,r be three real numbers <0 € and m € IN<set 

of +ve integers). 3 a +v© constant KC0,h,rJ> depending on O c 
T<0,h.,r> only through 0,h..r. and such that ^fCieTiO ,h. .r} -»■ 

is a linear and continuous extension operator, s.t. 

|P„| < KO.h.r) . 

We shall now state the Re 1 1 i ch-Kondrasov compactness theorem which 
has got important applications. 

Theorem 2.2 (Rel 1 ich-Kondrasov) : - 

jf 

Let O c R be a bounded open set of Class C . Then the 
following inclusions are compact. 

(i) if p<m., -♦ L^CCD ligKp* fwhere p* - 1 

^ p ^ ^ ^ m-p 

( i i > if p=n, >^<0> -► L^CCO li«?«» 

( i i i > if p>», -► CCCb . 
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The notion of trace has very important applications in the 
boundary value problems. The trace is a generalization of the 
concept of the restriction of a continuous function to a 
submanifold of its domain of definition. We shall now state the 
Trace theorem as follows. 

Theorem 2.? (The Trace Theorem for 

Let O be a smooth bounded open subset of D?’' and ao be its 


boundary. Then the trace operators 0 < h < m-i , can be 


extended to continuous 
Further the 
continuous linear mapping 


linear operators mapping onto 

operator y = (j'^.y^..., is a 

of onto n ^ Moreover 


a a continuous linear (right) inverse y ^ mapping 11? icsci> into 


The details about Sobolev Spaces can be found in the text books of 
Adams C19753, Oden and Reddy [19763 and Mason [19853. For the 
general concepts in functional analysis one may refer to Yosida 
[19743. Taylor [19583 or any standard text books in functional 
analys i s . 



CHAPTER III 


A SHAPE-OPTIMIZATION TECHNIQUE FOR THE CAPILLARY SURFACE PROBLEM 

5.1 : INTRODUCTION 

The study proposes a method to construct the capillary 
surface, based on the shape-optimization technique. The problem 
of determination of capillary surface is reformulated as an 
optimization problem over an admissible class of sets. An 
existence theorem for the optimal solution is proved and a 
numerical schem,^ is proposed based on the finite element method 
and the gradient method where the gradients are calculated with 
respect to the shape of the domain. Some model problems have been 
solved. 

5.2 : LITERATURE SURVEY AND GENERAL BACKGROUND 

Whenever a liquid and a gas remain adjacent to each 
other without mixing, the interface between them is termed as the 
capillary surface. Although the capillary phenomenon was observed 
even in the times of antiquity a consistent theory capable of 
scientific predictions was first proposed by Young C1805]. He 
introduced the notion of mean-curvature of the surface and could 
relate it to the surface tension. He furthermore observed the 
constancy of the angle of contact of a liquid surface with a 
solid. Laplace C1806] derived practically the same results as 
Young. The theory was further strengthened when Gauss C1850!l 
arrived at the same equations using the theory of virtual work. 
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The equation governing the capillary surface is the 
"Laplace-Young" equation which is a non-linear elliptic equation. 
The general solution becomes particularly difficult not only 
because the equation is non-linear but also the domain of 
definition is not known a priori. This is termed in literature as 
the "free-boundary" problem whose numerical solution is in general 
difficult to obtain. The basic problem in which the mathematical 
study is aimed at is the problem of existence <or non-existence) 
of the solution surface under various conditions. Perhaps the 
first published work of this nature is by Concus and Finn C19653. 
The paper discusses the existence of the surface in a domain with 
corners and the discontinuous dependence of the solution on the 
boundary which was a significant contribution. The 
discontinuous dependence has also been verified experimentally. 
Similar mathematical works on the existence were done by Concus 
and Finn C1974a.b.c, 19763, Finn [1974. 1979. 1983, 1986. 19883. 
Giaquinta C19743, Guisti [1978,19803. A comprehensive account of 
the mathematical theory is given in the monograph by Finn [19863. 
Though there are no generalized numerical schemes to solve 
capillary surface problem with volume constraints. in several 
practical situations some simplifications are possible. One such 
simplification is when the domain of definition is known such as 
when the liquid is contained in a cylindrical vessel and gravity 
acts along the generatrix. 

The construction of the capillary surface by numerical 
methods was studied by Concus [19683. and later by Siekmann, 
Scheideler and Tietze [1981a. b3. They used traditional methods 
such as expansion matching along with the finite-difference 
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methods. Since for fixed domains variational formulation is 
possible. several researchers used variational Finite Element 
Method which has the advantage of handling complex domain 
geometries easily. Such an approach has been used by Orr 
C1975.1977]. Brown [19793 and Mittelmann [1977] for the capillary 
surface problem. Brown [19793 used Galerkin Finite Element Method 
for solving the capillary surface problem of a liquid in a vessel 
of square cross-section. He compared the efficiencies of various 
finite-element. methods with different interpolation functions. 
Mittelmann [1977] derived some a priori error estimates for the 
Finite Element Method for the capillary problem. 

Another class of problems that are simple to solve are 
the ax i -symmetr i c problems [e.g. spherical and conical vessels 
etc]. In these case although the domain is not known a-priori, 
the dimensions of the problem can be reduced and hence the problem 
becomes simplified. A number of studies were carried out for this 
class of problems [Bashforth and Adams [18853. Chesters [19773. 
Padday and Pitt [1972,19753]. Although it is not possible to 
obtain a closed form solution when gravity is small, one can solve 
the equation explicitly for the zero-gravity case. The solution 
is expressed in the form of elliptic integrals and it can be 
represented in terms of Delaunay curves. However, use of the 
volume constraint poses some serious difficulties, as the volume 
is a quantity which has to be determined numerically and one does 
not know all the parameters in advance. Hence the explicit 
solution method cannot be used and it is a good choice to solve 
the boundary-value problem by means of shooting techniques. This 
technique was adopted by Bashforth and Adams [18853 which was 



possibly the first time when fairly comprehensive results were 
published and since then the work has become a bibliographical 
rarity. Ax i -symmetr i c profiles for ranges beyond the stability 
regions were computed by Padday and Pitt [1972. 1975}. A large 
amount of numerical data on ax i -symmetr i c equilibrium shapes is 
contained in the monograph by Hart land and Hartley [19763. 

According to the literature available the numerical 
solution of 5-Dimensional meniscus problem was reported for the 
first time by Orr et al . [19773. In their work Orr et al . used 
finite element method for the solution of the capillary surface 
equation. Their work was basically a simulation study without any 
attempt to justify the method. Under these circumstance there is 
a need to develop a numerical scheme with a proper mathematical 
basis and practical viability. 

The present work adopted a different strategy to solve 
the free-boundary problem. The potential energy of the system is 
evaluated as a function of domain and is minimized by changing the 
shape of the domain occupied by the liquid. Hence, it becomes a 
"shape-optimization" problem. For the numerical scheme the 
problem is solved using finite element method and a gradient based 
method. The gradient of the functional is calculated with respect 
to the shape of the domain defined by the nodes of the 

f ini te-el«fment mesh [Pironneau [198433. worthwhile to 

mention that the Augmented Skeleton Method <ASM) developed by 
Hornung and Mittelmann [1989a. 1989b. 19903 is fairly a generalized 
method to solve the capillary surface problem which takes into 
account the volume constrint and the surface need not be 

represent ible as the graph of a function. Though the A.S.M is a 



general method it is felt that the proposed method is more 
efficient for the class of problems in which the surfaces are 
r epr esent 1 b 1 e by the graph of a function. 

3.3 : NOTATIONS AND PRELIMINARIES : 

Let us consider a mechanical system comprising of 

liquid, gas and solid as depicted in Fig. 3.1. The liquid fills 

the vessel partially and the gas or the vapour of the liquid fills 

the remaining portion of the vessel. The walls of the vessel are 

considered to be rigid and sufficiently smooth. 

Let D denote the total region of the vessel O^the region 

occupied by the liquid and the region occupied by gas being 

denoted by O . The liquid has a specified volume vol . That is 

S 

the volume of D is finite and larger than the volume of the liquid 

in order to make certain that the liquid remains strictly inside 

the vessel. Let and Z represent the surfaces of the vessel in 

I € 

contact with the liquid and with the gas respectively. The 

interface between the liquid and the gas is denoted by S and y 
denotes the line of contact between the surface of the vessel and 
the surface of the liquid (Fig. 3.1) a point on which is denoted 
by P. 

Let denote the density of the liquid, the density of 

the gas is small and is neglected. The sur face-tens i ons of the 

liquid, liquid-solid and the solid-gas interfaces are denoted by 

O',, O' and O' respectively. 

I s g 

Let us denote (o' -o' )/o, by cose . G denotes the angle 

' Q S C 

of contact which is the dihedral angle formed by the liquid and 
the vessel on the line of contact. The liquid, the gas and the 
vessel material are assumed to be homogeneous. 
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The external forces acting on the liquid are only due to 
gravity and assumed to be independent of time. Thus the volume 
density of these forces is given by -pStn , where n is the 
gravitational potential. 

Thus, the total potential energy of the system can be 
expressed by the following 

O = "-ilSl + O-s |S:il+°'g l^gl + Pj ; " ■“l 

[where (.| denotes the measure of the corresponding surface]. 

Fig. 3.2 shows the representations of the system under 
consideration x^,x^ and « are the Cartesian Coordinates and g is 
the acceleration due to gravity directed as shown. The vessel lies 
entirely in the region a: > O. fi is the region formed by the 
projection S onto the plane 2=0, and P denotes the boundary of 
O. 

Let denote the region formed by the projection of the 
vessel surface onto the plane 2 = 0 , and P^ the boundary of Ci^. 
The surface 'S’ is expressed as a single-valued function fi) . 

The vessel is assumed to be open and its surface represented by a 
single-valued function yCjf) and yf € C*(Q^>. 

Hence we have for the system 

xi-yf>oyxeCi < 3 . 2 ) 

(from the conditio’ that the surface of the liquid lies entirely 
within the vessel). 


From the volume condition we have 
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FIG, 3. 2 notations 


J“ (u~yf)dO = vol <>.?) 

O 

Expressing the potential due to gravity as 

n = sz <5.4> 

the gravitational potential energy can be expressed as 

Ug = f n dOj < 5 . 5 ) 

Thus substituting <3. 4) in <3. 5) we get 

? 9 

Ug = -i — iu -yp -XKi <3.6> 

O 

Let Jk be the set of admissible domains O , which are 

•2 

non-empty closed and bounded subsets of (R defined as follows : 

J^. = tfi ; Q c O c O .. n .0.0 . satisfy the cone-condition 
a dad 

CTheorem 2.1 (Chapter II>]. The metric over is defined through 
the Hausdorff metric d|^< . , . > . For any Q e M let define the 

following : 

= <^l € <0) : C <C i s large and independent of O) ; 

w i y' V X € O and J dO = 

O 

Let be the unique solution in of a given problem over O 
and let G = C<0,u^>; O e . 

Let us denote ips/a^y by ho. 

Using equation (3-2). (3.3). (3-4). (3-5), (3.6) the energy 

expression (3.1) takes the following form' : 
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L/<a.-u> = J Cl + |v ^ dbc-cos 9^ f Cl + j'7 


(/)j ti+13 dx V 


■u e and O € ^ 


<3. 7) 


Let us denote L/<Q,u^) as J(0) 


?.4 PROBLEM DEFINITION AND MATHEMATICAL DETAILS : 

The equilibrium-surface problem can be treated as 
"opt imal -shape design" problem as follows : 

Determine Q € -41 such that 

-7<0*> < -7<Q) V O € ^ <P> 

"if* 

Then <0 . is called an optimal pair. For a general reference 

on the optimal shape design problems one may refer to Pironneau 
C19843 and Hasslinger [1908]. 

We now proceed to prove the existence of a solution of <P> with 
the help of the following theorems and lemmas. 

Theorem 5.1 : Let lO >6© a s»Qxionc& of a closodL and bounded sots 

ffi 

in Jk. Then 3 a ©ub«©<yu©nc© CO J> lubich. conx^oTg&s to a compact stot 

ft^J 

Q G M in tho HoMsdorff sonso. 

Proof : See [Pironneau (1984), p. 51 Theorem 1], 

Lemma 3.1 : Let I (u> = <l + |yu|^>^^^ + ~ Ju^ V ueW^ (O) . 

a O 

Let u — ► u in W^<Q) 


then 1 im / (u > = J<u>. 
n 


n-Mo 
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Proof : Let \> * 

„ n 


It is clear that v 


n 


-►V 


a 


2 . 1/2 


and I,<ti ) * 
2 n 


We have 1 im = I,<u> 

z n z 

rv #- ot > 


|/,(V„)- 7 <«)|= 


I V oVjrj 





^ J ll^nl' - l”l'l 

o 


===> lim /,<v ) = /,<v) 

1 ~n 1 „ 

n->» 

Thus the lemma is proved. 

Nom we make some remarks which will be useful subsequently. 
Let a cc n <Q satisfies the cone-condition). 

3 a linear map s <0> ► P^u = v <u « W^<0) s.t. 

P^j * u and also l^ol <K (AC is a constant independent of O) . 


The second assertion is due to Chenais C197?l. 

Lotnma 3.2 ; Let O ,0 cc O and w ,v € A/^ <0> s.t. v |_ e 

n n n ii 

n n 

Let V ► V in A/^(0) as Ci — ♦ fi. 

n n 

Then v|^ c V^. 

Proof : let us define the following 

( i ) w « 

(ii> wZ “ ■ <v„-v'))/2 

n ■ n • n 
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( 1 i i ) 


' T "n' 


t>ol ) ' 


O 


n 


<iv) = I j <w;;>^ 

O 

n 

Let D be a region defined as follows : 
m ® 

D € W and D cc Os.t.I> cQ Vn>n. 
m j m n m 

and D ► O as m. -► oo. 

m 

Clearly Q = D U C where C = D \i> , 
n m m m n m 

From Rellich’s Compactness Theorem we get 

s? 


n 

Thus 


V in <n) 


n 


-> in < 0 > 


-Xj) j s ^ I <v^-V/> - (v-v^> + 

- I 

^ I'—nl 




n 


-»■ w in (O) 


j, <w ) = i ~ vol)^ 

■'In 2 J n 

O 

n 

= i <J vol -J 

D C 

m m 

lim j <J >*> - - J 


rv+oo 


m m 

1 im lim <j'’j(u>^>)= y ~ w>l>^ = 0 

O 


n» 0 D n-»oo 


-v 2 


-'V'*n> ' 7 ‘f '"n’ 


n 


n 




m 


m 
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^ylim < >)=i f ® 

4L n c 

nvMso n-xx) ^ 

o 

="> "In ‘ ''d. 

L&mnai 3.3 : ThQ ^aph G carr^spondin^ to CFJ> is ss<pjtsnt tally 

compact in ths /olloxt>ing sonso : 

Let CO > <0 « M) i>s a soQuonco. Thsn3 a subssovtsnco 

n n 




and an slsmsnt COpU^J> € C s. t. 



Proof : First of all let us note the graph G is non-empty since 




<Pj> 


exists (see Baiocchi C19843. p.l3> in view of the fact that is 

1 

convex, closed and bounded subset of H Let the unique 

solution of <Pj) t>e denoted by (for convenience). From the 

definition of and the uniform extension property we get 

n 

S u I- -. < X. Thus 3 a weakly convergent subsequence Qu > which 

n." 1 ,ii H- 

n # 1 # % 

converges to \te (Q) . Let us define u^= Now we shall prove 

that is a solution of 

(P, > with O replaced by Q. 

1 n 
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let V ► V in iCi') with '^n . |_ e v_ 

’Pn. ^n. 

' ' J ) 


with the notations used in the lemma 2 we get from the minimality 

of in <Pj> 
n 


j [(' • I”;- iT" ■ 

n^ ' j ' 


> 0 


[(> ^ \T' -(' ^ l’“n lT'1 H- 

m ' ’ m. ' ' m 


( > > 0 


From the lower-semi cont i nu i ty and the continuity property (Lemma 
1 ) of the functional we get as n. ► <» 

I [{' ^ -(* ^ ^ I P - “o > ^ 0 


Passing to the limit as m -*■ a> we get 


i [(> M- -(• ^ " I 


Also as [from lemma 5.2Dw^ is a solution of (Pj>. 

~ 1 

L^mma 3.4 : Let Q — — — ► and ti u in // (Q> where v | € 

n n n|rt O 

|0 n 

» n 


and € V^. 


Then 


1/(0. ^t) < lim inf 1/(0 .u > 

n n 

n-M» 


Proof : With the help of the notations used in the previous 


lemma we have 
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= X[(‘ ^ ^ X ‘ ^ f it' 

m m m 

*(^]l (* " l!’“nlT''^ i. 


< . ) 


m 


lim <lim UiCi , v )) > 1 im 
rtH-oo n-^oo m+oo 


(a-'^T 

'■ m 


fl . feo 


+ ^ r <^ l -¥' ) 

m 


-cos 


®c I (‘ * x^'-’ } 

m d m 


===> lim uici , V y > uiCi.uy. 
n n 

n-K» 

Theorem 3.Z : Let the lemmas 5-3 and 3-4 be satisfied. The^ 3 
atleast one O , CO c which, is a solxLtion of CPJ>. 

Proof : See [Hass linger Cl 988], p. 29, Theorem 2.1]. 


3.3 THE DISCRETIZED PROBLEM : 


Let O be discretized to by triangulation is 

composed of 3-noded linear triangles CCiarlet C1978]] such that 


the vertices situated on F. also belong to F. We also have ,41. = 

n n 

CO. : O c O. c 0.>.ln such a triangulation we consider the 

n a. n d, 

h ^ h h h 

finite element spaces. and defined as follows ; 

o 

= CVj^ c W^<0^> n C*^<0^):v^ > y in O^ ; J Cv^-y) dx = vol;|v^|j< C 

^h 


'^h^r. € P,<T.) and v. |_ = Wr> s-t the nodes! 

j 1 j hlF^ F^ 


1/2 
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= £v. « W_<0. > n C <n. >: h'T. 1 J 

Cj^ h O h h J 

W^* = iw. « «l(n. > n C^CQ. >: '^hlr. ^ v/„ at the nodes} 

h O h h ) "h 


(where T e T. and P, <T.> is the space of linear polynomials on 
) h 1 J 

r. >. 

J 




h a *‘h 


and 


oh o 




J. <n. > - E. ■*• E., ^ A . 
h h h hi oh. 

Now the discrete problem is as follows : 

Determine such that < -7^(n^>V ^^h^ 

In order to determine Ej^ we solve the weak problem which is given 
below : 

Determine u. g v!? such that 


“h'^h * 


— <Pih<“h'''h” ■*■ 


Vu. 

h Do 


<Phl> 


'^u. . Vv. 

~ h h 


where a. (xi. ,v. >= f ^-tt? <*f + 6 o f xi. v. 

^ °h ” + l!“hl ^ "h 


dx 
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^h ^h 


dx) 




dbc 


and / =1 for ; and / = 0 otherwise ># = 0 for > y and 

W = 1 otherwise «g> 0 are penalty parameters. 

In the ensuing numerical scheme for the solution of 

require the gradient of -7|^<0^) . For this we develop the necessary 

tools . 

It can be seen that can be cast into a variational inequality 

formulation. COuvaut and Lions [19763, Glowinski [198133. 

L&ttma 3. 4 : L»t Tj tho *l»m&n.£ obtained, from. T j € by 

M M M M 

translating ons of its xtsr-ticss to <?’ ^ CFig. 

4%^ 0^ ^V* ^0 ^ 

5,3 Lot N*" C .!> bo tho basis function associatod. with, tho -oor-tox . 


A J* 


Thon Af‘ Cx^ - W CxJ> 


+ 0C\6q^[:> V x « T n T^. 
A# A# Ar J J 

In s 

V 9 ,<? vorticos of Tj 


<5. 8) 


Proof : Soo [Pironneau [19843. pg- i02']. 


Lomma 3.5 : Lot g bo a continuously difforontiablo function on T .. 

» h hk ^ 

Lot T . bo obtainod from T . by translating to q + 6q . Thon wo 
3 3 ^ 


hm>« 

T -CT nr :> 

J J J 




^ g,dx = ^6q^. '9Cghi'>* OC\6q*^\> <3. 9) 


r . 

J 


Proof : Soo [Pironneau [19843. pg. iOSJ. 
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FIG. 3.3 DEFORMATION OF 

T; BY SHIFTING o'f 
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Let us consider the variation of the solution when the node <? of 

T € T- shifted to *+* with the help of the following 

J ^ -V A. 

1 emnnas . 

Lomma. 3.6 ; Lot K donoto tho sot of hK>\mdjcxry ruodoo of and 

k 

tho ^sot of intorior nodoo of O,. Lot tho nodo q of T , bo shtftod 

^ -V 

k k 

to Q ^ Sq . Thon w hax>o tho followirtg : 


Cq/^y.6<f - 6q^ + OC |6<5i^p <?.10) 

Ai. A# ^ ^ - 


A^ A^ 


where X. is a characteristic function defined as follows 
h 

\ = 0 for 

A.. 

= 1 for q^G 

A# 

£ i. 

Proof : We have u. = E «./f + J) u .N 

^ ieK *" 


.6u. - r <<5u./f*'+ u. 6N^'> + r <«5u.W*'+ u. 6N^ ) 

h i i i i 

icK icAC 


<5. 11) 


From <5. 8) we get SN = 
(5.11) to get 




C<? J>.Sq which we substitute in 

A# A# A# 


t«AC i«ZiC -u teAC *■ n 

- E f^u/7t/cq^:>.6q^ <3. 12) 
teZC 


Ar A^ 


But Vti 


= E V <y .w") + E 

icAC * fcAC*" 


V <ti 


.//> = r u. + r u. vn’’ 

‘ " -V ~ 


<5. 13) 


Substituting (3.13) in <3.12) we get 


Su = £ 6x1. + f^X 7 V* (<?*>. 69 ^- Cq^:>. 6q^. 

icAC *■ -«.-v 
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&u 


We define sZ*! = £] //*■ = E 


t , h .,i 

1 ~ • ^ * 


Thus we get <5. 10). It is readily seen that c . 

L&mma. 3,7 : L^t r%odL& q af T . 6^ tron^lat^d to + 6^51 > iH&rt 

-A# 

thj& followtrx^ haldt V xi^ « and Vti^,6v 


,,h 

‘h*'^'^h * 


<i> ^ 


^7 C6u^:>.'7 

h h 

^ -V 


h 


dx - 


i 


V \i^.vc<5u.:> 

h h 

:r — cvu^. vv^;> dx 

h 


hO SxjL, xj.dbc 
n n 


+ r — — cvcyv^x^vy/p c«?^>.vv^- V 6<?^;>. Vv^:>dx 

•'• -V, h 


A# A(» A(» 


I- r — ^ vv^:> 

ft L A# A# Ar A# A» A^ A#- 

h ®h 

-cvv^.v c//^7^l^c<?*:>. 6<?^:> cvu^. Vv^:> I «ix 

A» A# A# A# Ar A# A^ J 


Vu. 


r _2l — . W, . 6<gi^J><ix + bo J‘//^CVti^C<3i^J>. 6 q ^ 


- Vti C< 3 «^^. cbc —bo jT dx 

A# A# A^ O A^ Ar Ar 


[[ 


CV^1- . Vv. ^ 
A A 

A# A^ 




j/Z'Jdx ^ 6<?\ dx 


+ OC j6<3>^p yihoro = 1 ••■ 

A» 


1 


C3. lA) 
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( i 


i> fiu* / 


/ dx) 


+ (/ < //'x.Vv'r C<f^:> .Scf^- //* 7 ri 6q^>fdx:> <r V / dx> 

ft ‘ h -w -V 

h n 

- i^jCvL^if)^/ dx - x>ol >C^ hi^ / V'Vj^C< 5 i^^. 6 <ji^cixJ> 


A* 


dx - vol )C^< 5 <y^. VC//'v^/J>cix^ 
h h 


VC/Z't 


+cf6<?^. V c/cwj^-¥»:>//^:>dx;K: f / v^«ix> + oc|6(?^|;> 

“h*" "h 


<3. 15) 


<iii) 


a<^2h<''h> •"h’ ' 1 " *“h "h'** * X'^'^h’^r. <<31* ) 

ft ft -V h 


-hi^xi^ic^y .6<^'>v^ dx 


^ -V 


I 


M <^l^-V'> f^v^iq^y .Sq^ dx 


-V *v 


+ J 6 q^ {M(.M^-yfyv^?/^ydx + 0 (| 6 «?*j 


O -V -V 

h 


Proof , <i) a^(u^.w^>|j.. = ^ — 


'^h -^h 


<3.16) 


dx 


I ^) + |Vu.| 


* ^ '‘h'’h 


dx 


V <\i.+6w. >- <v.+6u. ) 

n n h n 


<^K’K'>\Ti = J ?- 

h h h J) ^ (l+|S 7 uj^+ V 6 ri^|^> 




Air 


dx + r 6 o<rii.+ 6 'U^ ) <v^+ 6 v^ )dx 

J 




Vu. . Vv. 

h n 

Ao A# 

__ 


] dx + r 60 v^W^)dx + 0 <j< 5 <j^| 

2 * ^ ^ ^ 
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The last two integrals are obtained from <3.9> 


= J 

T . 


( V\i.+ Vv. ).(’7v.+ V6v. ) 
h h h h 

__ 


1+2 


Vu. . V<5v. 
h h 

-V -A# 




] 


- 1/2 

dx 


+ bo f V. dbc + 6o f u. 6v. cbc + bo f 6-u.v. 

^ h h f ^ 


dx 


i 


( W . . Vu ^) 
h h 


h 


1/2 


dx 


bo JT + 6q^.^ 


Vv. 
h h 


1/2 


] 


dx 


+ bo ^ (.u^'o^I^ydx 


Invoking (3.10) we get 


n h 


3 


3/2 


•-I 


r-T- V6u. >( Vtt. . Vv. >dx 


j 


* X -^2 <'»''> •■s?* 

/ . i . ^ A# K A- ^ 


i 


- V<A?N7v^<qi^) .6<?*>. Vv^ dx 

A#« A^ A» Aw 


^ -^2 [ K- } 




- f Vv^. <V^).6<?^<Vv Vv )>1 dx 

^ A# Ai» A# A# A<» Aw 
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X. *V --W -V -V -V ^ 


dbc 


( i i > 


-t>o .6qi^ <ix. 


n n 


+ f 69^.7 Y7T~ }}*^ X <tt^V|^A/^)<±>c 


But < 5 a|^(^lJ^ ,v^> = E • Hence we obtain ( 3 . 14 ) 

j j 


</ 3 j^(«h>-''h’lr. " 'I <'‘h‘»'’>' 

’ ^ j 

’ j 


<ix - VO I > < 


^ '’h ^ *'> 


<ix - vol ) ( 


I V' / dx) 

j 

dx - \>ol 


«(Pjh<''h>-’'h>lr’ = <1 <<''h-'">+ / 

^ j 

+ J Sq^ .'7ihl^<‘u^-yf}/>dx:>C^ (Vj^+ 6v^) / dx 

+ J 6<?^.V(/V^v^/)dx>- (v^-y»)/ dx - vol ) (^ v^ / dx) 

r. X. 

= <^6v^/ dx)<|v^ / dx>-(^/(v^-V')dx - vol ) ( p/'/ t^v^ <q^y .Sq^cbO . 


J 


+ <J^ /<u^-^)clx - XK>1)(^ Sq^.V (v^M^/)dx +(^ 6<ji^.V</Vj^-V/>A/*>dxJ> /v^dx) 


Invoking ( 3 . 10 ) and noting that 6 </?j ,v^> = £ .v^> 

j j 


We get < 3 . 13 ). 
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<iii) ' r 

^ ^ j 


* r V<«(u^-V)Vj^//')dx 

I . ^ -V* 






f 

1 


A## 

6 u* cbc 
n 


) 

i. A#- A^^ 

J 

-f M<u^-Y'>N\v^<<f^> .6<}^dx + r S^\vCM(u^-v;)v N^) dx 

jt . ^ ^ X . ^ ^ 


J 


As tV^) = 2 It <5.16) 

i ) 


Theorem 5.5 ; L«t w. .v* fc>© th» solution of P. , 

n A hi 

trinagulations and. T^ rospoctivaly, vihoro T^ is obtained 

JL^ jy 1^ 

by tr-anslating tho ‘oortox q to q' . Lot 

A# A^ 

solution of tho following problom. V Vj^ € ; - 


wi th 
from 
tho 
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i 


ITT- ^ T 7 ?— ^ X 

I. n. -%• -V O. 




h. e 




= -f cq^:>.6<}^:> - 7 c//Vu^c<?^:>.^*>). 7 wj^d 5 * 

jfj, '**^ ^ ^ ^ ^ ^ ^ 


•h ^h 


f -i-j—^CCVu^ 7 Ca/Vu^C<7^^.6<?*>)(7xi^.7<a/'x^7v'p 45<j*‘:>J>:>(7u^ .7v^)<Sx 


"h ^h 


Vu, 


7 u^. Vv. 
h h 


+ f -—n-i • 7C//* Vti. S<f^^dx - f Sq^. V CC 

A ^^ 1/2 ^ - - ft. - - 


•h ^h 


-bo JT A/^CX^T^p Cq^S>. Sq^J> - Tu^ Cq^J>. S<f^J>v^dx-*-bo f Wj^wN?Vj^C< 3 i^J>. 

Q ^ ^ n# A^ A» Q --W A-. A# 


dbc 


-ho r. 6 <y^. yCit^v^A/'^cix 

O. -V. 


( I /''h^ 

h h 

f r /C^l^-y:><^x-voll jjr a/' /7v^ C< 7^:>. 6<?^^dxj 

X fth 'ft^ -v -u -w 

^ ^ Cu^-V*:)/ <ix-volj ^\7 C//'v^/J>dxj 
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— f /ff/'cx. vy/- Cqf^:>.S<f^ - Vw C<?^J>. 6<7^;>v dx 

^B a . *^h -w 

n 


i-f XC^^ -y:>AA7v. 6<?^dx -^r 6<?^. 7C>#C•^^ -¥/J>N^v. :> dx 


(3.17) 


Then |U|;:^-ti^-6u* +//V C<?^:> . 6<?^-//'x^V Vj, c«?^;>. | j< C 


(3.18) 


wher« C is a constant vhich dspsnds on only. 

Proof : Combining (3.14), (3.15) and (3.16) we get (3.17). 
Using (3.10), (3.17) can be written as 




V6^t. . Vv. 

ix ft 

-<v -V 


1/2 


r ■ 

■ k- 57 ?— •" »° J^tTh'^ 


o . 








.VC//* Vu fij*J.<59*7<lJc 


1/2 
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-f <5o^. V^rr — -r-zrx fCu. v.C<y^S> . S<y^S>dx-ix:> f S<y^.VCzi v N^J>d 

L ^ M J ) L ^ ^ - 6 ..- - 


O -a* --u 

h 






/ dbc-'ool 




/ 7v. Cg^;) 




. 6<?^ dxl 




. . / Cxi-yp^ <ix-^>ol 

i A ^ 

h "h 




^ ' 7 C}^- 




^ L (^ <ix j 

h h 

r— f yC/fCi/ ~y/^^v -> <ix 

B w. ij ^ 


0 . 19 ) 


Let us define : fco = 5o + 


// 


’B 


and A^Cw, 

A 


h-'-h^ = 1 . 


7ii>, . Vv, 

A A 

-V -A# 




1/2 


CTw. . 7 v, J> C 7 w, . 7 v-^ 

A A A A 

_ Ar. -A. «A# "A,. 

dx - J. — TTZ~' ^ 


I 




bo w. V. dx 
n n 


^ i; (;[ -h / ^ ) [ X ''h / ) 


(3.20) 


It is clear that ./4. Cv>. , v. X> > O 

h h h . 


'= O • : ==>^ ^ :S=‘' ^ G ■ 

h h h n 
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1/2 H H 

hence qualifies as a norm. 

h 


Let us show the equivalence of j-jj and j • 


T T 2 

Define 7u. = Cv u J> and TJf'. = Cw -> for D. cc R 

h n, n„ n h, h. h 

^ 1 £ ^ 1 Z 


Substituting in <3. 20) we get 


V“h-“'h’ = X S/2 {“S-S/ “h, 

12 12 11 


+u. ii>. 

h^ J 


So dx * i-[ ^ / dx j [ ^ / dx ] 


f "~T 7 T“M 1 +<l+ti?'J> w? ~2u. tiu 1 ^ d. r *1: 2 

i lI ^i) hj h^ hj h^ hj h^ J <ix+r ba 

”h 


‘h ^h 


2 

^ ^ '"hjlj ^ <^1 Kl? 


Also we find the < 1 + (as ^ C) 






^7w. . Vx»>- 

tl A 


^ Kll 

n ey. 


From the equivalence of semi -norms we get 


<^5 Kli ^ KU ^ S l>“hlj 

n 


dx 



< c 




hh ^ ^ ^ ^ 


Theorem >.3 shows that a small change in the node induces a small 
change in the solution. 

Corollary 3.1 : L»t <5^* « «>• a salxition. of CJ.17^. Lot 


. 9 u . . 

© = E “i ^ ^ ^ 

-V ^ 


TV) donoto tho compononto of © 

* 


CwJWr# n. £« th® dtimorvoion, of cc R >• 

« 1/^ • 
h *^ 00 ’ 


Thon wo haroo t/v«« 


following: - V « 1^^* “ 


i 


V c©;^,:>.vv. 

hi h 

^ A# 

J72 


dx 


"4‘ 


^ Aw Amt 






dx 




h , — 


C<? 




9 u, 


h ^ 




A 


1 <<vu^.v c//* ^ -’‘h ^^--'-h':~h 

V2 .V I ^ " A, I 


Aw Aw 
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Vu .Vv 


cbc 


‘/i ^ h 


3 xi 


ha f N^CX^ CC}^^ - r<?^:>^v^ cbc * bo J 


K " 

^ A 


L 


&U 


h dx 


^ V- ^ X _5 

Cq J> dx 


LA 


^ ' }A:> dx - — 


-bo I - Ctj. V. ■ 
h h 




^ [ I /C.^-v.^<i.-vol] [ C,'':> 






C / Ctij^-v/;> dxj [j v^ / dxj 


— r J MfAcX. C<}^:> - CqA^v.dx + — J MCu -^f^:>jA ^ CqA 

ejj I “v ^^7 J ^^7 '^o “ h dx, ^ 

B ^ h I 5 I , 




i = 1.2 n 


<3. 21) 


S h I 


Proof : This is a straight forward application of 


<5^ = E 6^i N* = E —A N^.Sq^ into <5.17) 


■mJ 0 ; 


■dx 
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Theorem >.4 : L»t 6 ® a solution of tho following probl 0 f»- 




Ai» *\i» 


* I jl* Pk / ‘=*^ 

h 






Vu. .Vv. 
h h 

■A# ■-%» 

TTz 


A ®h 




cix + f 6 o u.v.cbc V v. c k' 


h h h Oo 


( 5 . 22 ) 


># 

with bo = 6o + — 


TA®n + 2 ^;^ I = 1,2 n 


9 J’ 


9qt 


k. ‘hi h 2 h 3 


< 5 . 2 ?) 


where = f < Vu^. g f - Vu ^ dx 

n, '*•'»' l'»' -V-*. l'»' 


'h ^h 


^ * '^lA 


bo r 9 ^^^ 2 


^i":>cix 

n 


^ ^ dx 


krz 


9 u^ 


+ bo f<ti^ ^ jgJ^ y dx 

^h 


^^'^h " ^h " " '*h " 


7 ^-‘ = - CO, 


0 ^ ^ <aAc 1 + | 7 v/|^>^^^) dx 
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r *?7 ’"h 

O jT £ A-- -V t ^ ^ 


’h ^h 


t 


Vz <^'^h- ^ f <j''^:>)-<vu^. V g f J^cvu^. Vp^:. 

^ ->v-^ L ^ -v*a*> I ^ 

h ®^^ *' 




dp 


'^h- ’'•h 


-r ~i/2 • VC//* jgj^ fg*:):) dx -r ^ <c:;: ^^;^n'';> dx 

Ci. g~ I I g^ 


*h ^h 


n. l I Q^ I 


■*^ L 3^''h'^v ^ *J7'/ '*') lj5^'’hH 


" ^ i “fcj V *• ) 


' ^ U C'‘h-i^'^<'>‘ ) ( I '>h / H 


dx 
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i— J Mh^CX ^ ) p dbc 

h 


+_ 1 _ r MCu.-^:> <q^> dx ~ — r -^MCxJL,-y>> pj> 


dx 


Proof 




dx 


-cos 9 


^ Cl + j^|^>^'^^dx + 


Oh 


.‘^J <0^>= c- .6<^ = f -2! ‘^—K — dx 


■h' 

«9 


k 


h ^h 


TTT 


k 


^^. 6u. dx 
n n 


^ 6c}^. ^Ci^e^^:> dx + ^ ^ <5<?^ . (VCW^t/^J> - VC//'v'^0)dx 


-cos 


et Sq^. vc/^ci+|^/v'|^->^'^^ ^ 

O. -v -V 'W 


h ^h' 

==> ■ V S<i^ * ^ — T7z — ^ * X ^^*^h ^ X 


9q 


h ^h 


Vu 


^ f - VCW^y/' 

fit. -V -V --V 

n 




i^chi^x^^ c<f^:>.6<f^ 
g . '*' n 


h '^h 
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S<f^y ]dx 

^ ^ ^ 


4 - 5 o r 




//'cx. Ctf^^.Sq^ - ^u,Cq^:>.SQ^:> dbc 


h h 


h '*' 


-V -V 


-cos ©^J <S<?^. WW^Cl+|^|^:>^'^^)cbc 

jQ A# A# -\» 


9 J 




i 


'’l “h "h 


172— “h '*=' 


i . 


9 u 

^ c ,*:. - ^ 

^ -w I ^ -W I 


h ^h 


•r fc .-»-. r > *. . !? X 


I - d * 


*bo dx - CO . e {• c ,^ cx .\^\ h^'h 

u ^ I ^ I ^ u I 


k 


^ u ..^ ce ?,^ 

n nt 

A# 


A ^h 


1/2 


i 


-h ®M ^ 2 ‘ 


Substituting in place of in (3.22) we get 


^ ^ f TT- ^ ”h ®h£ *' 


dbc 
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^ 4 [V n# ^ ^ « 

®h 


"X *>» <l Ph * 5^ li 

“h h 


Now we invoke Corollary >.l to replace the R.H.S. of <?. 

also substituting p. in place of v, in the Corollary >.l 

fX rt 

0.25) 


/ (X ph / J 

h 


3.6 DESCRIPTION OF THE NUMERICAL SCHEME : - 

Let ij - 1,2 no. of boundary nodes) denote the 

describing the position of the interior nodes of w. 

boundary nodes . 


i . e . 


k. r I.,, w A 

Q = G<C<?>>V<? 


e and e K. 

A# 


Now we can obtain the expression for the derivative of the 
functional as follows : - 







1 


9J^ 


«? 

Ar 


q c 

A# 




Oq 


ij * 1,2. 


.n> 


The algorithim is as follows 

(1) Choose C initial mesh], nmax, eptol 

n 

do( . >m = 1 . nmax 

(2) Compute from 

(3) Compute from <3. 22) 


<3. 24) 

24) and 
we get 


function 
r , t the 


energy 

(3.25) 
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9J^ . 

<4) Compute the gradient - - — j- ; <?*"« AC : i = 1,Z 
from <3.2>>. * 


,n 


(m> 


(5) Compute the approximation of of arg min J, ^ ^ 3 


V g « K 


4.1 w o 4. -(m +1 > . <m+l)l <m>l . (m) <m> 

<6) update the mesh Compute O. by g =g +p g 


w ^ ^ mf 1 y ^ > X 

V g € AC ;g = & Cfg >J> 


4 ^ 


n 


<:7.> if 


J Cg^^^^K »ptol stop 


t = 1 


C. j> Continue 


3.7 MODEL PROBLEM : - 

We consider a simple problem of the determination of the 
equilibrium-surface of a liquid in a conical-vessel. 

The following properties of the system are used. 



= 1.0 

a 

c 

o 

o 

It 

°'s 

= 0.001 

H 

= 100 

a 

= 0.8661 

vol 

= 100 


g 

a = 45° 

[where a is the angle of the cone and M the height [Fig. 3.4a]3 

3.7.1 Ax i -symmetric Formulation s 

The problem is an axi-symmetr ic one and hence it is 
reducible to a 1-dimensional problem involving the radius as the 




• • 
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variable. The finite-element discretization is performed with the 
help of one dimensional linear elements. In this case the mesh is 
kept uniform in which the coordinates of an internal node is 
directly proportional to its node number and hence the function 
relating the boundary nodes to the interior nodes turns out to be 
simple. The figure with finite-element discretization is shown 
in Fig ? . 4a 

The following Z cases have been studied for the ax i -symmetr i c 
problem. 

<i) 60 = 0. 

The abscissa of the contact point can be determined using the 

analytical formula given below CMyskis C19873]. 


R = 
c 


'• n cot OK 


il - 


cot oi 


sin <oi-d ) 
c 


(l-cos(oi-© >)^(2+cos(ot-© 
c c 


= 4.92 

This agrees with the results obtained using shape-optimization 
technique . 

The Z-coordinates of the computed surface is given in Table 3.1. 
The value of the gradient of the energy functional has been 
determined from (3.23) at various points and it matches with the 
one determined by the difference scheme. 

< i i ) ho = 0.1 

The Z-coordinates of the computed surface is given in Table 3.2. 

In this case also the gradient of the energy functional matches 
with the one determined by the difference method. _ 

The solution surfaces for the two cases is shown in Fig. 3.4b. 



O.OOOEOO 


0.426E01 


0.492E00 

0.427E01 

0.984E00 

0.429E01 

0.147E01 

0.452E01 

0. 196E01 

0.376E01 

0.246E01 

0.443E01 

0.295E01 

0.430E01 

0.344E01 

0.458E01 

0.596E01 

0.468E01 

0.442E01 

0.479E01 


0.492E01 


0.492E01 




Table -3.2 


Node No. 

R-absc i Bsa 

Z-ordi nate 

1 

O.OOOEOO 

0.432E01 

:2 

0.488E00 

0.433E01 

3 

0.976E00 

0.434E01 

4 

0.146E01 

0.436E01 

5 

0.193E01 

0.440E01 

6 

0.244E01 

0.444E01 

7 

0.292E01 

0.450E01 

8 

0.341E01 

0.457E01 

9 

0.390E01 

0.466E01 

10 

0.439E01 

0.476E01 

11 

1 - - 

0.488E01 

0.488E01 


* I - ■ r ^3 


ifiK. Ho. 


J,]402:i, 

I * f an* ' 
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1 n 
the 
The 


3.7.2 Two-Dimensional Formulation 

We now solve the problem described in (i) 
2-Dimensions . From the consideration of symmetry 

computational domain consists of a quadrant of a circle, 
initjal mesh is depicted in Fig. 3.3. There are 31 nodes and 43 
elements. The final solution surface is depicted in Fig. 3. 6. In 

this case the function relating the boundary nodes with the 
internal nodes is more complicated. In the present case the 
axi-symmertic property o the problem is exploited and the nodes 
are constrained to move along the radial directions. This 
considerably simplifies the data structure to be developed. The 
limitation on the node movements is put only from the 
consideration that the triangles should not be flat. The 
conjugate-gradient method has been used for the problem. 


3.8 


DISCUSSIONS 


An alternative method of the solution of the equilibrium 
surface problem has been obtained on the basis of 
'shape-optimization' theory. The major advantage of the method is 
that an accurate and inexpensive technique has been developed to 
compute the gradient of the non-linear functional. Hence it has 
become possible to use an efficient optimization algorithm based 
on the gradient for the numerical solution. In our proposed 
scheme the convergence is obtained in the sense of the 
optimization algorithms. It is to be noted that the technique for 
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the gradient computation is fairly general and can be applied to 
other areas of mechanics involving constrained variational 
functionals. As an example we may cite the 'Constrained-obstacle' 
problem which is of special interest to the researchers. 
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FIG. 3.5 2-D FINITE ELEMENT MESH OF THE 

LIQUID DOMAIN. 
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CHAPTER IV 


A FINITE ELEMENT ADAPTIVE STRATEGY FOR THE CAPILLARY SURFACE 

PROBLEM 

4.1 INTRODUCTION 

This chapter describes an h-version adaptive strategy 
f or ' the finite-element scheme and the computer simulation of the 
capillary surface problem. An error indicator criterion based on 
the residual error estimate is derived and verified numerically. 
A suitable alternative error indicator has been proposed for 
adaptive purposes. A model problem is solved using the two 
different error indicators. 

4^.2 BACKGROUND AND LITERATURE SURVEY : 

The shape of the capillary surface is governed by the 
"Laplace-Young” equation. When the domain O is fixed the problem 
is equivalent to that solution of a non-linear partial 
differential equation in two dimensions. Some amount of work has 
been done on the computational aspect of the problem. In the early 
stages finite-difference methods have been used by many 
researchers CConcus C19683 .Si ekmann C1981 a.bll. 

As the capillary surface equation is the Euler equation 
corresponding to the energy functional, one can apply the 
variational finite element method for its solution. The finite 
element method has been used by Brown [19793. Orr [1975. 19773 and 
Mittelmann [19773- In many cases the solution contains sharp 
gradients and may be inaccurate in those regions. So there is a 
necessity for improving the accuracy of the solutions. In the 
present work an adaptive scheme is proposed to solve this problem. 
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The general objective of the adaptive refinement scheme is to 
automatically change the mesh , or the approximating structure of 
the computational methods so as to improve the quality of the 
solution. Implicit in this goal is the resolution of two basic 
issues: the first one is to obtain a criterion for the assessment 
of the quality of the solution and the second one is adapting 
another mesh to improve quality. 

The first question is resolved by obtaining an estimate 
of the error of the solution in some appropriate norm. Since the 
exact solution is in general not known a priori an attempt has 
been made to estimate it by deriving a-posteriori error estimates. 
The error estimation based on the estimates of the residuals is 
used widely in finite-element analysis. The residual error 
estimates were introduced by Babuska and Rheinboldt C1978 a,b] for 
linear elliptic problems and have been applied by many authors. 
But for the non-linear case it cannot always be guaranteed that 
the error is bounded by the residual. That is the reduction of 
the residual in some given norm does not necessarily imply the 
reduction in the error of the solution. However, for certain 
classes of non-linear operators the error is bounded by the 
residual. For monotone operators it is well known that the error 
is bound by the residual COden C1986a3]. In this Chapter it is 
shown that the differential operator of the problem is strongly 
monotone and hence the error is bounded. 

As the estimate of the error is global in nature, the 
question arises as to how it can be used per se as a local 
refinement criteria. In general, the global estimates have been 
used as a basis of local enrichment of solutions. For the 
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refinement process to be efficient there is a need of a reliable 
local error indicator. It has been found that the error indicator 
based on the residual error estimate does not function 
satisfactorily due to the nature of the operator associated with 
the problem. But one can derive some other estimates based on the 
residual error which can perform as better refinement index which 
will depend on the nature of operators one deals with. In the 
present case a new error indicator has been derived which 
performed in a better manner as a local refinement criterion. 

The second issue is related to different adaptive 
schemes which are a subject of much research. For finite-element 
approximations three broad categories of adaptive schemes are in 
use: In 'p methods' the local degree p of the polynomial of the 
element shape-functions is increased to obtain a better 
approximation CBabuska Cl 981 a.bll. The most common method is the 
'/v-method' in which the mesh size h is reduced by successively 
subdiving an initial mesh thereby reducing the solution error 
CBabuska C1978 a.b], Oden C1986]]. The third type is the moving 
mesh method in which a given mesh with a fixed number of elements 
is fixed and the nodes are shifted so that the local errors are 
reduced. In the present problem the /i-method is adopted. 

Four-noded rectangular elements have been used for the 
numerical study. The data structure has been developed in which 
the original element is subdivided into 4 new elements COden C1986 
bll. As the governing equation is a non-linear one the 
finite-element formulation leads to a set of non-linear algebraic 
equations. The 'full' Newton’s method with block-storage scheme 
for the numerical solution is adopted. 
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4.5 PRELIMINARIES 

2 

Let O cc R be an open set and T its boundary which 
is Lipschitz. 

Let I' be a space of functions defined as follows 
y = M^<Q) n (Q> 

Let us represent the elevation of the surface by the function 
u<x.y) and let 0^ be the angle of contact Csee Fig. 4.1]. 

In this system the equations governing the capillary surface 
are the following 



bo 


9 


in O 


<4.1) 


V u . n 

— = cos 9 on r 

yl+jvu| ^ 


<4. 2) 


where n denotes the outer normal to Q and q, bo are known 
constants, feo i 0 . For zero gravity bo = 0. 

The variational formulation for <4.1> and <4. 2) can be written as 
follows CMyskis C1987]] 

Find u c V s . t . 


< r,v > 


- I 


O 




•/l + jVuj 


+ bo 


I 


cosB V = 0 < 4 . 5 > 

c 


Vv « k 
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We now consider the discretization of (4.3) using finite-elements 
Let be a triangulation composed of Q' (bilinear 

quadrilateral) elements AC, X « cr in such a way that all the 
vertices situated on the boundary r. of the set O, = UK also 

belong to the boundary. With such a triangulation let us define 
the space V’ <z V as follows 


V: = C ti. € W^(0. ) n (O, ) .* ri. 1 € Q' (X)> 

AAA A 

2 o 

From the imbedding theorem we know that — ► C (Q) 

Hence the interpolation is well defined. 


Thus for the discretized problem we have : 


Find s . t . 




J 




— h — A 


-/Tk 


4“ ho u. V, 
A A 






-I 


COS& V- = 0 

c K 




(4.4) 


4.4 A-POSTERIORI ERROR ESTIMATES 

Let us define the following function spaces 

(i ) = C vj « I * ^ ^ 

I ^ 

where P iK) is the space of polynomials of order p > 1 on K. 
P 

(ii> = C vf « V'f and K' interpolant of = 0 > 

ho A A A A 
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We have | r^|| ^ = sup > 

*' |u| < 1 

An approximation of the residual is constructed according to the 
following : - 


l“o ' ^ 

I'-fl ^ > 


where C is a constant, v € K and . 

o h h. 


If h. stands for the mesh size we have in general |u^— ^ 0<h.^> 
so that it makes sense asymptotically < as A — ► 0 ) to 
approx imate 

sup <r^ , V > by sup < r^, > 

Let us define the error - u [where u is the true 

solution] We have from (4.4) 






~ A 




7 u 


/ l + |7u^| ^ / l + |Vuj^J 


. V ©, + bo o, 

^ n, A 


Let V Vu = T 7 


Define <|> ,r>) = 


i / 1+K I 


/l + |r?P, 


<?-») 


K-r)| 


We find that <|>(0,7>> = 


> 0 


[where | . | 
stands for 
the Eucl i dean 
norm] 


/ l + hl ^ 

We now proceed to prove that <}> (? .t)> > 0 for ^ ^ 
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[/T+j^ £ - ^ ■/ ] ■ <f-^> 

4> <?.r?) = " " 

/T+| 7 i^ /i+\v\^ \l-v\^ 

' (KI^ /T+j^ + \r)\^ / 1 + |?|^ ] - (•/ l + + /I+i^ ] 

■/i+KI^ / 1+hl^ 

[l?l^ / i+h|^ln| Vi+|?|^ ]-Klhi[A7|^+A7]^ ] 

■/ -/ 1 ?-^ 1 ^ 

[|fj-|^|)(|?| /i+hp-hl/i+Kl" )(lll +|^iVi + |?|^ ) 

/i+Kl^ Vi+hl^ |?-^|^[l?ivl+i^|^ + 

< Kl + hi > < Kl -hi 

/i+|?|^ /i+hl^ K-^|^[|?|/i+hl^ + h|-/i+|?|^] 

The equality is valid only when ? = W for some « R . 

The only case to be seen is when % - ~ xi . 

* |f|^ v' 1 + |?|^ 

In that case <k(^ ,ri> = > 0 

< . ) 

Now we consider the limiting case when ? — ► >> - 

«w ** 

Let ■ T7j + p CO* 0 


and 


^2 * ^2 ^ ^ 
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<|> <P.©> 


2 1/2 

/ — ; — r 2p<77,cos©+T7p»ine p ') ■! 

yi + |i7| •{ (T 7 .co«d+T 7 ®tne)-<i7 co«©+r>-sin0) 1+ 5-= + 5- >■ 

~ I ^ ^ I 1 + hl^ l + |r)|^J J 


6 > / i + i?r ^ 1+hl 


=> <|> (p.©) 


1 + ? 


( 

l + lr)r ^ 


(77jCO£?9 + 

1 * i»i^ 


+ €<p>| 


where €(p> — ► 0 as p — ► 0 


Xl4h. < <|><p.©)) 
p-*^0 


1 


Tf^cosO + rt^sinO 

1 - _J £ > 0 

1 t |v,|2 J 


Hence we obtain for t>o > 0 

2 2 2 

< r, > > a|«, }, _ + bo le. In _ 2: c, I©. I, ^ C for some constant 

ri n. ' n.' 1 ,Q- ' rt'O.n. 1 * n"l,n, 

'* '* c, > 0 ] 

'> Kl-i.o, Kl ^ '1 Kli.n, 

. A n, 

/V 

When bo =: 0 we consider the quotient space V = V /y <y is a constant). 
We have < > £ “I'hli.o^- 

Also from the equivalence of l-lj q I'll O obtain 

* H * h 

h h 


Thus the residual gives some indication of the error in the 


global sense. 
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To estimate the residual we solve the auxilliary problem ; 

Find w, € s . t . 

fx no 



where is the solution from the neighbouring element. 

^ (U I! “hi']’'' 

where C = max C. [ C. > 0 ] . 

^ h h 

2 

The local error indicator is = J* i.e. we can use 

K 

o. as an estimate of the local error over each element. In 
h 

general reducing implies a reduction in q which in 

h 

turn implies a reduction in ju - j q • Although the 

' h 

residual norm is the global indication of error it can still be 
considered to be a criterion for the refinement. 

4.5 THE PROPOSED ERROR INDICATOR : 

The simplest strategy for adaptive refinement is to 
compare the error indicators of each element and refine the 
elements which have large errors. In general the regions of sharp 
gradients will have large errors manifested in the form of large 
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error indicators. But in the type of problems considered here the 
effect of the gradient is not pronounced in the error indicator. 
This is due to the fact 


V u . n 

that in the boundary integral we have whereas in 

U + |V u\^ 

ordinary linear elliptic case it is Vu . n . It is clearly seen 

that the modulus of the error in the boundary integral tends to 
get normalized for large values of hence the error 

indicator of the above kind does not satisfactorily function as an 
index of local refinement. 

In view of this we propose a new error indicator which indicates 
the regions of sharp gradients in a better manner. 


Let 


fru 


denote the jumps across ^K. 


We define the following 


< '■ch • ''h > = 


V . 


Z'‘k 


■li + 12 


+ t> U. + Q 
oh ^ 




< ‘oh • ^9Krf. 



1 


sr 


- cos 0^ 

.1 

•1 J + l^^hl^ 


&Krr, 
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It is possible to compute the constants Cj^ 


and C 


and hence 
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ermirted. However^ due to the computational 

“ 1 , fC 

difficulties we use the following simpler error estimate suitable 
for practical purposes. 

We again define an auxilliary problem 


Fl/id <j>^ « <X) s.t. 


k -tic Z t * i < ff ^71 ■ < ^ 


« x \ r . 


Therefore we have 


V vf € 

fX fX 


< 


■ ‘ Ik ix-z "h 


X X 


= x> 


ii'^Aii-1 ^ = (1/ 1 n 


1/2 


We notice that 


fr^n 


0 as h — ► 0 


= = > 


9n 




0 as ^ 0, 


Therefore 


^ (2-1 I 1^ ) 

A 


as h —»■ 0 . 

1/2 


as h — ► 0 . 


We now define our error indicator as e. 


= J I ’tc I' • 
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The practical advantage of using e^ is clearly demonstrated in 
the model problem. The corner regions near the vessel wall have 
sharp gradients for the solution and hence maximum refinement is 
done in those regions CFig. 4.5b]. 

4.fr ADAPTIVE MESH STRATEGY : 

The important steps in h.-ref i nement strategy are as follows : 

Step 1) Form the triangulation out of O. 

Step 2> Solve equation <4. 4) using non-linear iteration scheme 
e.g. Newton’s method. 

Step 3) Compute the error indicator for each element K 

and form a refinement index [where (j is 

the mean and & the standard deviation and 1 > ^9 > 

0.1 

Step 4) Obtain the solution at the new nodes using 
interpolation and the Lagrangian multiplier approach to 
take into account the constraint nodes. 

Step 5) Repeat from step 2 until ‘s small or the maximum 

number of elements permissible has been reached. 

4.6.1 Details of Data Structure in the Computer Code : 

An important consideration in all adaptive schemes is 
the data structure and the associated algorithm needed to handle 
the changing total number of elements, nodal locations, numbers 
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(b) 


(c) 

FIG. 4. 2 SEQUENCE OF REFINEMENTS 
OF A UNIFORM MESH . 
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and element labels. The data structure developed in the present 
problem is similar to the one developed by Oden et al . C1986 b] . 

We shall now briefly describe the data structure which has been 
developed: 

In the present refinement process an element is refined 
into 4 sub-elements. The book keeping of the elements and the 
refinement process is monitored by the arrays NELCON <8, NED. 
NELSID <8. NED. XCOOR<NUMNP) YCOOR ( NUMNP > . and LEVEL <NED 
where 

I'HJMNP is the number of nodal points in the system. 

NEL is the number elements. 

XCOOR is the array containing the X coordinates and 

YCOOR is the array containing the Y coordinates of the nodes 

NELCON is the element connectivity. It is to be noted that 

although we 4 noded-quadri laterals are usedthere can be 
4 additional nodes on the sides due to refinement. Thus 
we have to keep the information for 8 nodes in an 
elements. 

NELSID is the array containing the information of the 

adjacent element. 

LEVEL is the array which assigns a level number to an element. 

Initially, the same level is assigned to all the 
elements. For example the initial coarse mesh could be 
assigned level 1. When an element is refined the 
sub-elements belong to a higher level (2) result. This 
process can go on indefinitely in general. 

The algorithm for refinement is described in the following steps. 
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St«p t Loop over the neighbours of the element < which is made 

possible with the NELSIO array) and check the level of 
refinement of the neighbouring elements ralative to the 
level of the element to be refined. 

Stop 2 If any neighbouring element has a level lower than the 

element to be refined then refinement is not posible at 
this state. 

Stop 3 If the element can be refined, we generate new element 
numbers and modify the necessary parameters. 

Stop 4 Compute the modified NELSID and NELCON arrays. 

Stop 5 Adapt the connectivity of the negihbouring elements 
Stop 6 Interpolate the solution at the new nodes. 

It is clear that some strategy is needed to test if the 
elements are appropriately connected for the refinement to take 
place. Let us consider, for example, the uniform grid of four 
elements shown in Fig. 4. 2(a) and suppose that the error estimater 
dictates the element A is to be refined. Thus, A is divided into 
four elements. 1. 1 1. 1 1 1. IV. as shown, and the solution values at 
the junction nodes, (shown circled in the figure), are constrained 
to coincide with the averaged values between those marked X. It 
may also be noted that the connectivities change in this process, 
e.g.. the connectivities 4 and 8 of element B are different. 

Next let us assume that an additional refinement is 
required, and that we must next refine element III, We impose the 
restriction that each element side cannot have more than two 
elements connected to it. Thus, before III can be refined. 
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element B must first be refined, as indicated in Fig. 4.2(b>. The 
constrained node P in Fig. 4.2(a> now becomes active, while node 
Q remains a constrained node. Mith B bisected, we proceed to 
refine III into subelements cu.ft ,y ,6 , and new constrained nodes, 
again circled in Fig. 4.2(c), are produced. In this case, only 
elenient B had to be refined first in order to refine III, but, in 
general, the number of elements that must be refined in order to 
refine a particular element cannot be specified. 

4.7 NUMERICAL RESULTS 

4.7.1 Numerical Experiment 

We perform a numerical experiment to verify the 
assertion that the norms of the residual and that of the error 
decreases with refinement. For this purpose we take up a simple 
1-D problem whose solution is known. The problem is as follows :■ 
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As usual the V'^ <0^) space is composed of linear elements. 


The space is composed of quadratic polynomials (i.e. p = 

Z). The norm of the residual is computed according to the 
def ini t ion 



sup < V, > 


A constrained optimization algorithm is used for this purpose. 
The variation of the norms with refinement is given in Table 4.1. 
It can be seen that the residual norm gradual decreases with an 
increase in the number of elements which vindicates our assertion. 


4.7.2 Z-Dimensional Example 

We now apply our adaptive scheme using the two error indicators 
to solve a 2-Dimensional problem as shown in Fig. 4.3. 

The problem is to determine the capillary surface over a 
rectangular domain in zero gravity. 9^ = 45° makes the 

gradient infinite which makes the solution difficult. The 
solution is arrived by two methods : 

I) Refinement with error Indicator : 

The initial mesh < 10X10) is shown in Fig. 4.4a and the final mesh 
after 5 refinements are depicted in Fig. 4.4b. 

The variation of the error indicator norms with the number of 
elements is given in Table 4.2. It may be observed that gradient 
is high in the top right corner and indicator norm the error 
decreased by a factor of 10 in 5 iterations. 

ID Refinement with Error Indicator . 
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Table 4.1 


No . of 
Elements 


Error 

Indicator 

Norm 


4 

. 300 

E-05 

.255 

E-01 

.257 E-01 

9 

.755 

E-04 

.155 

E-01 

.118 E-01 

14 : 

. 556 

E-04 

.729 

E-02 

.792 E-02 

19 

.189 

E-04 

.557 

E-02 

.594 E-02 

24 

.121 

E-04 

.425 

E-02 

.475 E-02 

29 

.840 

E-05 

.515 

E-02 

.396 E-02 

49 

.502 

E-05 

.208 

E-02 

.257 E-02 


Table 4.2 


Refinement No. of 

Elements 

Error 

Indicator 

Norm 

1 

100 

4.250 E-02 

2 

150 

5.944 E-02 

5 

196 

3.559 E-02 

4 

400 

2.154 E-02 

5 

487 

2.012 E-02 







FIG. 4,3 COMPUTATIONAL DOMAIN. 
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The initial mesh (2x2) is shown in Fig. 4.5a and the final mesh 
after 5 refinements is depicted in Fig. 4.5b. The variation of 
the error indicator norms with the number of elements is given 
in Table 4.5. It may be observed that the elements are densely 
clustered near the top right corner which shows the effectivity of 
the, error indicator. 

4.8 DISCUSSIONS : 

An a-posteriori error estimate has been developed for 
the non-linear problem which is also verified numerically. An 
alternative error indicator has been proposed whose practical 
advantage is also shown. The h-refinement procedure seems to be 
effective as the time spent in the logic of refinement is less and 
the solution is interpolated at the new nodes which ensures fast 
convergence at higher level. 
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FIG. 4.4a INITIAL MESH 
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Table 4.5 


Ref inement 

No . of 

E 1 ements 

Error 

Indicator 

Norm 

1 

4 

0.228 

E-00 

2 

15 

0.144 

E-00 

5 

25 

0.981 

E-01 

4 

67 

0.551 

E-01 

5 

82 

0.514 

E-01 

6 

105 

0.485 

E-01 

7 

155 

0.448 

E-01 

8 

265 

0.277 

E-01 

9 

516 

0.258 

E-01 

10 

588 

0.245 

E-01 
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FIG. 4.5b MESH AFTER REFINEMENT WITH 



CHAPTER V 


OSCILLATIONS OF AN IDEAL LIQUID IN LOW GRAVITY 


5.1 INTRODUCTION : 


In this Chapter a description of a numerical scheme using 
finite elements is presented for the computation of the 
eigen-values of oscillations of an ideal liquid under Low gravity 
conditions. The major difficulty in this study is the fact that 
the differential operator defined on the boundary is of higher 
order than that of the spatial operator in the domain. 
Nevertheless the problem is cast into a standard eigen-value 
problem with a suitable choice of operators. A Finite Element 
procedure has been described in detail and numerical results are 
presented. 


5.2 LITERATURE SURVEY AND GENERAL BACK-GROUND : 

Linear oscillations of an ideal capillary liquid 
represents a classical problem in fluid mechanics. Lamb C19453 
has obtained some results in this study. The study of the liquid 
oscillations is of special interest in space-technology. The 
satellite contains liquid propellants and the magnitude of the 
natural frequencies of the propellant and its behaviour is of 
prime importance for the stability investigations. 

The dynamic movement of the liquid surface is termed in 
the scientific literature as "sloshing". A substantial amount of 
study has been carried out on the sloshing behaviour of the 
liquids under normal conditions Budiansky C19603, Bauer C1965a.b]. 
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Levin C1963]. Tong [1966], Zienkiewicz et al . [1978], Haroun 

[19813. But not much work was done about sloshing in 
micro-gravity or zero-gravity conditions. Bauer [1981, 1984, 1989 
a,b,c] adopted a semi -analyt i cal method for the computation of the 
eigen-frequencies. In the work vessels with simple geometries 
were' considered and this method becomes very difficult for the 
complex geometries or even in simple 3-D problems. 

A comprehensive treatment of the lumerical methods for 
the computations of the e i gen-f equenc i es is given by Myskis 
[19873. But the numerical methods covered do not include the 
Finite-Element Method. An outline for the construction of 
operators so as to reduce the problem to a standard eigen-value 
problem was discussed. In the present work this concept has been 
extended to the construction of the operators in the Finite 
Element spaces with a mathematical justification. Also a proof 
is furnished on the existence of the eigen-values with 
appropriately chosen function spaces. 

9.3 NOTATIONS AND PRELIMINARIES : 

A system comprising of a liquid, a solid and a gas is 
depicted in Fig. 5.1. The liquid is assumed to be the inviscid 
and i ncompr ess i b 1 e . The region occupied by the liquid is denoted by 
Ci. The surface of contact between the liquid and the walls of the 
vessel (assumed to be rigid and smooth) is denoted by E 
equilibrium surface is denoted by r. SP denotes the line of 
contact of the liquid surface with the vessel and 9^ the angle of 



V V \ 
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solid / 


Liquid with 
Density pj and 
Surface tension 




FIG. 5.1 


Notations. 
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contact of the liquid surface with the vessel. 
We define the following function spaces. 


Vj«| < j. 


wh^re a is a multi-index, and m > O 


is any positive integer. 


vcco = j u = oj. 

o 

vcr:> = |v€W^cr.>; J ^ = oj- 

r 

= trace of vt on r 

For conveni ence^, in further discussion we shall use ^^p in place of 

For the details regarding the space id) and the trace operator 
one may refer to Oden C19763 or Adams C19763, 

Let denote the velocity potential of the fluid and X the 

eigen-value of the system and X = X 

Governing equations in the fluid are given by : - 


=0 i n O 

3d> 

^ - O on E < no- flow condition) 


} 


(5.1) 


where n denotes the outward normal to O on £ 
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The following Euler equations govern the f ree-sur face . 

-•V ^ ^ ‘X’ ^ ^ '‘r r . 

with the b.c. 

9 ,9dt^ 94> ‘ <5.2> 

~ 9(k •> 

and the volume condition J ^ = O 

r 

Where t>oC^^ for ^ 'for s € 9r are given 

functions depending on the system properties. c is a constant 
corresponding to the Lagrangian multiplier for the volume 
constraint and ® is the tangent vector Csee Fig 5-1] 

We assume 6o € CD and x « L^C9rj> with x 

occuring in (5.2) is the Laplace Beltrami Operator on the 2-D 
manifold P. (the details are given in [Appendix 13). 

The distinguishing feature of the boundary value problem (5.1) 
(5.2) is that the equations featuring in (5.2) contain derivatives 
of higher order than the spatial equation in (5.1) : 

The present structure is not readily amenable to the variational 
method of solution, but can be reduced to a standard eigen-value 
problem with a suitable choice of operators. We shall proceed now 
to define the operators. 

Let u denote the normal velocity of a point on F. 


Thus ^ 
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p p 

Let us define an operator B : L cr;> — ► L CrJ> as follows 

p 

For a given € L Cr> l»t ^ denote the solution of the following 
problem 

Find 4> € VCO> s. t. 

J Vv = J w <5.5> 

Q r 

V V € VCQ> 

Then Bu 

The properties of B are given in the Lmnuna 5 .1 

L&nma 5.i : - B is a compact , sslf -adjoint and positive ds/inits 
operator-. 

Proof We note that the norms ^ and ^ are equivalent 

for [Deny C195533. Also the trace operator -y^i (0>— ►L^<r > 

3 

is compact for OcdR CCiarlet C198833. Thus B is compact. Let 

2 

be a solution of <5. 5) for a given (D and a solution of 

2 

<5. 5) for a given veL <r>. 

Then we have 

JCBu>v= JV Hence B is self-adjoint. 

r Q r 

Also J <FuJ>ti = > O : Hence B is positive definite. 

r n - - _ 

We define an operator M : H^CF^ — ► L^CF^ as follows : — 

C/#u,v:> = CSu.wJ) = I £CBu,%>:> + CBv,%03 



We define an operator A 


H^cvy 


h^ct:> 


as follows 


CA-a,x>:> = J7^C^l,v:> 

r 


Jbou. V + J* x*^v VveW^ CP^ 

r ar 


Here is the Beltrami first differential operator on P 

CMyskis C1987], Nitsche C1975]]. 

p 

When ho > O clearly for ti€KCPJ> CAu,xO > ot P for some a > O 

<5. 4) 

When bo<0 inequality <4) may not hold in general . We restrict 

ourself to the case in which <4> holds CHI. 

2 

We define an operator G: L Crj> — *■ VCF^ as follows : 

2 

For any g c L CPJ> , Gg^CT^ denotes the solution of the following 
problem : 

Find Ge^cr:> s. t 

CACGg:>,M:> = cg,-oy ( 5 . 5 ) 

V ^^(P) 


2 2 

Lemma 5.2 :G is compact as an opsrator- from. L CPJ> — >■ L <rp;> 

Proof From Re 1 1 i ch-Kondrasov Compactness theorem H^(P) is 

2 

compactly embedded in L CP.>. Hence <? is compact. 

2 

We next define an operator T = GoB : L CPJ> — ► V'CPJ) as follows: 

2 

For a given ^^€L CPJ) let Tu be the solution of the following 
problem : 

Find Tu e /Cr.> s . t. 

CACTu:>»v:> = CBu,v:> (5.6) 


Vv€ KCP.> 
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p p 

Lomma. 5.3 : T is compact as an operator from. LTCT^ to L CrJ>- T is 
also sol f -adjoint and positive dsfinits. 

Proof B is compact from to L^crj>, C is compact from 

to z.^<r>. 

Hence T is compact 

Defining <xi,Tv>^ = CBxi, T\>> we have 


<u,Tx>>^ = J^^CTxi,TvJ> + JboCTvL>CTx>^ + Jx^TxOCT-oJ> = <v,T%i>^. Hence 

r r ar 

T is self-adjoint. It is clear that CBm,TxO > O. Hence T is 
positive definite. 

From the theory of compact, self-adjoint operators CTaylor [1958], 
Yosida Cl 9743] we find that T has a countable number of positive 


eigen-values with 0 as its cluster point. 

Let « be the eigen-value of T and W the e i gen-vector . 

m m 

TW = u W 
m. '^m. m. 


= T Cp ^ W :> 
m m m. 


-i 


^ W = T CX W :> [where X = P i 
m mm m m 

Thus X is an eigen-value of the following problem 
Find C‘u.x:> « vcr:> X fR* s. t. 


CAutVJ> = X CBu,v^ 
= X CWu,v.> 
Vv^cr.? 


<5.7> 


We can summarize the above results in the following theorem. 



91 


Th&or&m. 5.1 : - 

The system given, by C12> and. C2:> has a countable number of positive 
eigen-values xthich can be arranged in a non-decreasing sequence 
i .e. < O < < X^. . . and an orthonormal basis of 

eigen- functions w.r.t the inner product 

5.4 FINITE ELEMENT FORMULATION : 

For OcdR'^ let 7^ denote the triangulation of O into using 3-D 
finite elements. 

Similarly for TcdR^ let denote the triangulation of F into F^ 

using 2-D finite elements. 

We define the following function spaces. 

h 

^ •• 1 “h = ° ^ 

h 

where P CKD denotes the polynomial of order m. restricted to K. 

fh. 

We wish to make some remarks on the modelling of the surface. 
Since F is a 2-D manifold 3 a di f feomorphism G which maps and 
element AC on F to a standard finite element AC in the plane 

CFig.5.23. Let P denote a generic point on F. Then P has the 
following representation. 




FIG. 5 . 2 Surface representation . 
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= Cx^C?.r>;>. 


, >>:> , x^cf , T7;> 



£,<pt N^Ci^,r)^ denote the shape-function associated with the node i 
of an element and n the total number of nodes of the element. 


n 

Then x. = £ . 7 >:>x . . J = i ,2,3 

J i=i J'- 

Ulhere x.. is the value of x . at node i. 

J 


We also define the following tensor components. 



2 



S 

; P = 

. r 

9 ■^p 

X d - p 



9r) 


2 


G 


* r 
9 -^p 

9n 




CE6-F^>^^^ 


The Beltrami first differential operator is given by 


5.4.1 Construction of Operators and : ~- 

Let be the solution of the the discretized problem of (5.5). ie 
Given 8.t. 
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J ? ? = S '5 ®’ 

Using standard finite element interpolations for and we 

obtain the following discretized equations of (5.8) in matrix form. 


' e sa 

where denotes the array of nodes of and the array of 

nodes of Let denote the nodes of on 

Thus we have 


«r ' tr e 


= a »L 

s: — ^ 


(5.9) 


\j)h»r» = Zp ^ ~ ^p- 5 

Here H is termed as the coupling matrix. 

We have the following weak form of (5.2) as the eigen-value 
problem 

Find « R* X *■ t 


J * 

f-A 




* I ^'V'h 


= X 


A J 

r. 






(5.10) 


cr.i 

A A 
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<5.10) can be expressed in the matrix form as the following : 

ir 

A 

" = «•« 

= !ih>'h^2ih Where ^ 

Hence the discretized eigen-value problem corresponding to <5. 7) is 
the following : 

Find C\^pU^S> « X R^ s. t 

^ =?v -h < 5 . 1 1 > 


Mhere N is the resultant dimension of the system. 

(5.11) is a standard eigen-value problem which is also symmetric 
and positive definite and standard softwares are available for its 
solution. The operators ^ and ^ correspond to the stiffness 
and the Added-Mass operators under normal conditions. 


5.5 NUMERICAL ALGORITHM 

The major task in the numerical solution is the 
construction of the discretized operators 4^and We shall now 

describe some details of the implementation of the method of their 
construction : 

St«p i The first task is to construct the free-surface. If it 
is not known a priori, the problem of the construction 
of the capillary surface was discussed in details in the 
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Chapters III and IV. Once the surface is constructed we 
perform the discretization of the domain and the 
f ree-surface . The surface need not be discretized 
separately. The domain is discretized using 

3-dimensional <20-noded) solid elements and the sides of 
the element which lie on the surface of the liquid form 
the 2-dimensional discretization of the surface. The 
information regarding the coupling of the fluid nodes 
and the f ree-surface nodes are required for the 
construction of Tj, matrix in <5. 9). 

St»p 2 We now construct the ^ matrix given in (5.9) by solving 

(5.8). Instead of directly inverting stiffness matrix in 

T 

<5. 8), we can store the 4 § ^ factorization and solve 
for a load matrix which consists of the identity matrix 
to obtain the inverse. The advantage of this procedure 
is that some standard finite element method softwares 
can be used for static solution in a compacted storage 
form. The matrix H is stored in a secondary storage 
unit, 

St«p 3 The next step involves the computation of the matrices 
A and B. . They are computed at the element level. The 
matrix is similar to the ordinary stiffness- matrix 
except the fact that some terms from the boundary also 
contribute to it. The matrices are computed on the 
basis of eq” (5,8). It is to be noted that fully 

populated and not diagonal which occur in the case of 
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lumped mass matrix. 

Stop 4 After the matrices and are constructed, we 

directly invoke some eigen-solvers for positive definite 
matrices with compact storage mode. Since many 

finite-element softwares available for the dynamic run 

use the lumped mass matrix approach some significant 
changes are required to be made in the module to account 
for the fully populated mass matrix. 

5.6 NUMERICAL RESULTS z 

The proposed method is now applied to compute the 

fundamental eigen-values of the free-surface oscillations of an 
ideal liquid in a cylindrical container of square cross-section, 
under zero gravity conditions for various angles of contact 
(Fig. 5.?>. From symmetry considerations only a quadrant is used 
for the computational purpose. The 3-D mesh is depicted in 

Fig. 5.4. 

For the domain occupied by the fluid 20—Nodod solid olomont 

CZienkiewicz L197711 and for the fluid surface 8—Hodod 

Qxiotdri latoi'cil are used. A 3—D mesh generation routine has been 

developed for this purpose. The computed values of the fundamental 

eigen-values for various angles of contact is shown graphically in 

Fig 5.5. It may be noted that for the above geometry the 

free-surface does not exist for & < 45° [Finn [198631. For this 

c 

minimum contact angle the eigen-value turns out to be negative. 
For higher contact angles positive eigen-values are obtained. The 
eigen-values show an increasing trend with an increase of the 
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System parameters : 

Density of the liquid r 1000 

Surface Tension ' = 72.75; xio^ Kg/s'^ 

(gas -liquid interface) 

FIG. 5. 3 Computational domain details. 
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FIG. 5.4 3”D Finite Element Mesh. 
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angles of contact up to 90°, This trend is similar to the results 
published in CBauer C1989c]] for the liquid oscillations in 
rectangular containers. 

The results can be compared with the theoretical 
eigen-values of the oscillations of an ideal liquid in a 
cylindrical container of similar dimensions and with 90° contact 
angle . 

The fundamental eigen-value <X) is 0.7E-05 units for the 
latter system. The order of magnitude of the eigen-values are 
almost same for the two systems which enhances the authenticity of 
the results. 

9.7 DISCUSSIONS 

With the construction of operators ^ and ^ thereby 
reducing the problem to a standard one a numerical method using 
the Finite Element Method has been successfully applied to solve 
the eigen-value problem. As one has to solve two different 
problems for the construction of operators the use of parallel 
processing becomes attractive. With this method the dimension of 
the problem is also greatly reduced. This is due to the operator 
B.. One can also study the fluid-structure interaction by adding 
the structural stiffness and the mass operators with little 


difficulty. 



CHAPTER VI 


OSCILLATIONS OF LIQUIDS IN VESSELS WITH TUBE BUNDLES USING 

HOMOGEN I Z AT I DN TECHNIQUE 

6.1 INTRODUCTION : 

. In this chapter a method to analyse the oscillations 

of the an ideal liquid taking into the effect of capillary 
surface by applying the homogenization technique is presented. 
This is a study of the asymptotic behaviour of the eigen-values 
and eigen-vectors in a periodically perforated domain with the 
number of holes increasing to infinity. The distinguishing 
feature of this problem is that it is a 3-dimensinal one with a 
periodicity in Z-dimens i one . The method of asymptotic expansion is 
used to obtain an explicit formula for the homogenized 
coefficients for the problem. Using the method of energy 
theorems on the estimates of the eigen value and their convergence 
are proved. Numerical examples are presented. 

6.2 BACKGROUND AND LITERATURE SURVEY : 

In the last few years some efforts have been made to 
study the vibrations of tube-bundles immersed in a fluid. Such 
problems arise in the mechanical behaviour of heat-exchangers, 
condensers, reactors etc. This problem has been numerically 
studied by Chen C1977], Paidoussis tl9773. Schumann [1981], 
Shinohara [1981], Tanaka [19811. Planchard [1980,1982,19851, 
Planchard and Ibnou Zahir [19853 and many others. 
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In practice one is required to solve this kind of 


eigen-value problems where the number of tubes is very large [e.g. 
forty thousand in the core of a P.W.RD. If one uses finite 
element method, it is necessary to have large numbers of mesh 
points to represent the fluid between the tubes leading to 
matrices of large order (eg 10^ for 2-D problems). The cases 
become worse for 3-D problems requiring still larger memory which 
may not be available on general computing machines and large 
C.P.U. time. Thus, there is a need to identify a simplified 
problem which is equivalent to the original problem. One such 
simplification is the limit-spectral problem which is a result of 
the asymptotic behaviour of the original problem. One such 
technique is the Homogenization technique. 

Homogenization deals with the study of the limiting 
behaviour of the differential operators of physics in 
heterogeneous or perforated materials with a periodic structure of 
periodicity € as « 0. For a heterogeneous material the material 
coefficients are not constant but vary with the space coordinates. 
In the case of a material with periodic structure the coefficients 
vary periodically. If the length of the period is very small one 
may think that the solution of the differential equation is 
approximately the same as the corresponding solution of a 
homogenized material with constant coefficients. The 
homogenization method, shows the existence of solution and derives 
some properties of the homogenized operator. 

The theory of homogenization in perforated domains has 
been developed by many authors. The mathematical aspects of the 
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problem have been studied by Lions C1980, 19811. Duvaut C19763. 
Cioranecu and Saint Jean Paulin C19793. The eigen-value problem 
in perforated domains has been studied by Vanninathan C1978 a.b, 
19813. It is mainly the results of the works of Vanninathan which 
are useful to construct the homogenized operators in the problem 
of vibration of tube-bundles. Planchard C19803, Planchard and 
Ibnou-Zahir C19853 have studied the problem of vibration of 
tube-bundles using the theory of homogenization. 

6.3 PROBLEM DESCRIPTION ; 

Mhile considering the problem of the oscillations of 
the capillary surface of an ideal incompressible liquid the 

following assumptions are made : 

1) The walls of the vessel are rigid and enclose a finite 
region in R^. 

2) The tubes are rigid and are placed in a periodic manner. 

?) The motion of the vessel takes place in low-gravity 
env i ronment . 

The present problem of tube-bundle vibrations is 

different from the problem studied by the earlier authors. The 
earlier studies neglected the effect of f ree-surface , which is 
justified under terrestrial conditions. Also the tubes were 
attached to flexible supports which provided the potential energy 
to the system while the liquid was considered to be 

incompressible. The problem also could be reduced to a 2-0 one 

from the considerations of symmetry. In the environment of 
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low-gravity the effect of surface tension becomes pronounced and 
the study of the oscillations of the capillary surface becomes 
important. As the frequencies of oscillations of the capillary 
surface are much lower than the structural frequencies, the 
structure is considered to be rigid in the analysis. The elastic 
energy of the system is provided by the capillary surface which 
behaves like a stretched membrane and the kinetic energy is 
provided by the liquid motions. Thus the problem is a essentially 
J-D one. This problem can be thought of as an extension of the 
study of the oscillations of the capillary surface in an ordinary 
vessel . 


6.4 NOTATIONS AND PRELIMINARIES 

Let us consider an open set' 6odR^ defined as follows : 

G = C(x^,X 2 ,x^>€<o, > x<0, l 2 >x<o, l^)} 

T 5 

(x^.x^.x^) is a representation of a point in R . 

dimensions of G along the 3-di rect ions . G contains a large number 
of identical rigid tubes placed periodically in the x^ ,x^ 
directions. The axes of the tubes are parallel to x^ CFig. 6.13 
The boundary of G excluding the top surface is denoted 
by and O = represents the top surface of the domain G. 

represents the walls of the vessel. 

We define the following set M 

/f ® £ (x^ ,X 2 .x^)«G; x^« .constant} 
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where M represents the cross section of domain <Fig 6.2). 
sake of convenience we define x = 


For the 


Let 7 




T 


a 

and ( a f 

Let y denote a representative cell CFig 6.3] and T the hole. Any 
cross-section (normal to x^> of the domain is composed of the 
translation of the cell Y along Xj and x^, directions. 

We define the set Y* = Y/T 

-Let €>o denote a parameter associated with the periodicity of the 
system. The periodicity of the system will be « and « along 

the x^ and x^ directions respectively. This is also termed as c Y 
periodicity. 


We define 9 


meas I Y* 
roeas \y 


The characteristic function X<y} of Y is defined as following 

X CyJ) = f V y €y* 

= O V y €T n y 
2 

X<y) is extended on 0? by periodicity. 

The holes T*<js=l,2..> are defined as the closed connected 
components of the set 
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o 

O 

0 

o 

o 

o 

O 

0 

0 

o 

O 

0 

0 

o 

o 


dJI 


FIG. 6 2 A VIEW OF THE CROSS-SECTION 
OF THE DOMAIN. 
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FIG. 6. 3 DESCRIPTION OF A 

REPRESENTATIVE CELL V. 
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Cx; XCxy^s^ - O. X e O} 


N 

„ y € 

We shall denote M = M/ U 

* / j =l ^ J 


where denotes the number of tubes in the model. 

N 

' * / * 

Similarly n =Ci /u .T® 

« /j=f ^ J 
N 

and G = G / U .T® 

« / j^i ^ J 


where in this case T^. represents a tube correponding to 
T*. We shall make the assumption that the holes do not 
dCi. We denote by my.</) the average of / over /. 


myC/;> 


and my^Cf:> 



|y| y 

f /Cy;>dy 

■'y* 


Let us define the following function spaces : 


vca^:> = {v€W^cq*.> ; f v = o] 

q" 

vca> = iy>»H^CCD ; / V = O } 

O 

k<r6*:> = cv«-w^cG*:> ; f v = o } 


the hole 
intersect 
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ycG :> = cww^cg:> ; J v = o } 

G 


li^= Cv®w^cy''j> : V is y periodic} 


6.5 MATHEMATICAL FORMULATION : 

4b £ 

Let ^ denote the velocity potential of the fluid and 
the velocity of the capillary surface movement in the Xj 

direction. Let n denote the outer normal to the domain under 

consideration. Since the walls of the vessel and the tubes are 
rigid the normal component of the fluid velocity on the walls and 
the tube surface is 0. 

We denote the density of liquid by p, the 

surface-tension by o and acceleration due to gravity by g. Let 

bo = pg/a. We assume bo ^ 0 with bo = 0 when g = 0 (i.e., under 

zero gravity conditions) Let denote the eigen-value of the 
system. In the liquid domain we have the following equations : 


« _ - 
4* ' = O 

a 


trt G 


9 4> 
9n 


O art 9& 


'A* 


€/ € 


(P-l> 




on. 



o 
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Along with the existence condition = 


and the uniqueness condition J" = O 


On, the free-surface the following equations hold. 


-A^^ 


feou* + c = 




= o 


in O 


on Xi 


(P-Z) 


Mhere c is a constant corresponding to the Lagrangian multiplier 
for the volume constraint. Since the admissible functions satisfy 
the constraint condition we can eliminate the constant c from the 
equat ion. 


6.^.1 Variational Formulation : 

In order to transform the problem <P-l>+(P-2> into a 
standard eigen-value problem in the weak form we construct the 
following operators. 

We define ^ ^ as follows. For a given 

let 4^ denote the solution of the following weak problem 


Find K<G*> s.t 

J Y * = J 


V v€ vcg":> 
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Then ® ^ 

Ci 

c 

Me next define an operator CO* J> -> CO* J> as follows : 

C € 

<i4n^,t>^ = J" Vv + J* bo u^v Vvt*,v € CO*^ 

o- o 

We define A^: as follows .* 

O' <B 

= Cb“,v:>= 4 lCB^u^,vP-KB%,u^:>} 

With the help of -4^. and we obtain the following eigen-value 
problem of <P-1) + <P-2) in the standard weak form. 


Find C^e®.X®:> € VCO^:> X fR* s. t 
vJ> = vJ> 


=X^<A^u^.v> V v« 


CP-3) 


4E 

From the results of chapter V we note that O < ^ ^ X^ ^ 

where X^ is tho eigen-value of <P-3>. 

We intend to study the following problem. 


Find CxI^,\^^^CO:>X R%. t 
CJ4^^^^,v.> = X^C5^u^»v.> 

= V v^cn> ^ 

0- 
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where A ,E> and are the homogenized operators obtained as an 
asymptotic limit as « > 0 of the operators and <P-4) is 
termed as the homogenized problem of (P-3> 


6.6 ASYMPTOTIC EXPANSION : 

^ We now apply the method of asymptotic expansion to 

determine the homogenized coefficients for the problem. The 
method of asymptotic expansion was applied for the first time for 
homogenization problems in two dimensions for the perforated 
domains by Lions C1980. 1981]. We adopt a similar method for the 
3-D eigen-value problem. 

The asymptotic expansions of ^ and X* are the following ; 



We have A = € + 2 € 

Also we have ^ ^ lj-i»23 

Here A^ denotes the Laplacian w.r.t x variable and is the 
normal at point 

From (P-2) we have 

-Ati^ + feoti* = X*^^ i n 


(6.4) 
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Now substituting and X* from (6.1), (6.2) and <6. 5) in 

(6.4) and equating the 0 or -ve powers of c to 0 we get 


-A u = O 

y «> 

-A M.-2L. u = O 
y i xy o 


-A xi~2h u.-A XJL + b u = X d> 
y 2 xy i X o o o tx o 


From the boundary conditions on T we get 


9u 

o 

9n 

-y 


= O 


9u. 9x1 


9u 9u 

C *r-D = ° 

-y J 


From <6. 5) and <6. 8) we obtain 
ri Cx,yJ> = w CxJ> 

O O “ 

From <6) we obtain 


A XI. - O 

y i 


in y 


We introduce a function X*^<y> as the following 


-A = O 

y 


- « 

■ ar " — • w\ , 

9n j 

y 


in y 


<6.5> 

< 6 . 6 ) 

<6.7> 


< 6 . 8 > 

(6.9> 

( 6 . 10 ) 


(P-5> 


on T 
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X-^ iB y periodic, j‘ * 1,2 and n^. the component of the 

outward normals to T. It is easy to verify that 






< 6 . 11 ) 


Frorn. <6. 10) and (6.11) we obtain 


^ xk 9u. 

ifn “ ^ L ax. 

-y *■ j I 


on T 


< 6 . 12 ) 


Thus from <6. 7) and (6.10) we obtain 


6 u -A ti_ + 2-^ 
o o y 2 ^ . 


ax. ax . 
i J 


1- 7 u 

•' X o 


Vo 


<6-15) 


From the existence condition of solution of <6.13) we have 


au 

/ ■ J 

y' r y 


txJ r*'‘o , ...» 


J ^ Isrfc:] -I’' I Vo - l^' IV. 

V* i *-7 


<6. 15a) 


Since J nXyi:> = O we have from <6.12) the following 
T 

9u 


/ ^ o ; [ srSrl 


r y 


I J 


.th 


Let us denote by v a vector whichhas the i component unity and 
other components O. 

Then from <P-5) we obtain 


J 7 X^.v = J w n 
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Cwhere we have used the relation .v=7. <X-^v>- V.v and 


then the divergence theorem while noticing the fact that V.v 
since v is constant] 

Therefore 


= O 


r 9 X^ 

' «y. 


J 

T 


r 

J 'Sir 

T -y 


r 9xJ’ r 1 

J 3^ I ~^x :9x J 

_ L i jJ 


9 u, 


Substituting J* by J* 


9X-' 


r 1 

i 


<6.15a> we get 


9 X'^ 9 ^u 


r* o o X o 


-i— ( / )[ 

|y'| / '• 


0 

9x 9x . 

1 J 


] = 


X 4p 
cro 


Hence we have the homogenized equation 


9 M 

~Q • . ^ 

J 9x .9x . o o 

^ I J 


X ^ in a 
o o 


<6.14) 


where <? 


^3 


L— r 7 CX^-y J>. 7 CX^-y.y 
V*l J ~ *■ 


<6 . 14a> 


\y \ 


are the homogenized coefficients. In obtaining we have 

used the fact that 


J 7 <X'^-y - J X^VCX-^-yj3.n = O since 7 CX-^-y^.n = O on T 




If we denote by ,■3,. ,._i 7 


as the matrix of the homogenized 
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coefficients, then we find that C<7.,3 is positive definite. This 
is easily seen as follows. 

The condition for Cqi. .3 to be positive definite is 


"’^22 8ymmetric3 


But <?j 2 


-i— J Y <Ar^-y,)-Y <x^-yp < ( -^- J Y 

|y*| y* ^ ^ jy*| y* ^ 


< 


\y*\ 



<x^-y 2 ^\ 


2jl/2 


Cfrom Cauchy-Schwartz inecquality3 

We now try to express in terms of so that (6.13) reduces to 
a standard eigen-value problem. 

Using the asymptotic expansion <6. 2) in <P-1) restricted to we 

obtain 


-2. .< 
€ A o 

y 


+ 2c 


-i 


A A +A ^ + 




= a in y 


-Ay^^Cx.y.x^ = O 


-A^^^Cx.y.x^ -^A^^^Cx.y.x^ -A^^Cx.x^ 



O 


(6.15) 

(6.16) 


From the boundary conditions we obtain 

^p-<x.5j.x^ + n 35^x,x^ = O 


on r 
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d4>^ 

'SfT' 

-y 


9u 




on T 


as 


9x. 


= O 


on Q 


Using X'^C-^J>C j=i ,2> from (P-5) we obtain from (6.15) 


as 


j ^ xrJ^ % /''*'<> Cx,x^. + ^.Cx,x^ J = i,2 

S^Cx,y,x^ = - X-'Cy:> Cg^ - 3 •) '*'1 - 3r 


From (6.16) we obtain 


~\^2 ■ ^^xy^i “ ■ 




9x, 


- O 


(6.17) 


From the existence condition of (6.17) we obtain 


9X^ 


2 J IS^ C 9x.9x .] 


\y^\^o-s ^ 

T -y 


3^4> 

9xt 


A<p 


afr^fer^) 

f a “> 1 


\y I V* i J 




9x1 


= 0 


Hence we have the homogenized problem corresponding to (P-l> 


-o . . o 

a^r^. 

^ J 


a^s a^i>^ 


- O in G 


9x1 


94, 

9n 


= O on E 


^ “ ““o 


(P-6) 


on O 
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Where coefficients are the homogenized coefficients defined in 

<6.14) 

We define the homogenized operator in the following way 
Let be a solution of <P-6) 

^olo = 

Hence we obtain the standard eigen-value problem 






9x . 9x . 
I J 


+ m. „ (bo) 





t.J = 1. 2 


<P-7) 


where the coefficients . are given in (6.14a) 

6.7 HOMOGENIZATION THEOREMS : 

Let us define the extension operators P and Q as follows : 

P€ SiCN^cy*^:>, H^cY:>y s. t 

where C > O is a constant. & is the space of Bounded Linear 

Transformations. From Cioranescu [19791 we know that such an 

extension exists. 

" 2 » 2 

Cte ft CL CCi^^pL Cn^) is defined by extending the function 
to O value in the holes. 

It is to be noted that the problem <P-4) is similar to the Neumann 
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eigen-value problem, the only modification is due to the presence 
of operator ^ in place of the identity operator. For this we 
study the asymptotic property of in the following lemma. 

L»mma 6. i 

Let Pxj^ — ► tt* in L^CCD 

We have > in 


^ ^ ^ 

where a u = ^ I - where ^ is obtained from the solution of the 

O 

following problem 


E <?ij 

i)=1.2 

3 34>* 

1 1 

.2 * 

+ 1,^ = 

in 

G 



€ 

on 




•^3 

on 

O 


where q>. . are the coefficients given in <6.14a> which are obtained 
by the asymptotic expansion. 

Proof : Since Pu is uniformly bounded, from the definition of B 

we find that is bounded in L^(0> , hence~we can extract a 

weakly convergent subsequence. We have to show below that the weak 

« 

limit turns out to be 9^ . 

Define ^ .= i = i,2,3. In view of the uniform boundedness 

of jare uniformly bounded and hence 3 a subsequence s.t 
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We consider the weak form of (P-1) as follows 
Let DCO . Then we have 


G G ^ ^ a 

e € e 


+ J 0? . 


3 ax. 


= j 


V 


a 


Passing to the limit we obtain 


J i*- *J f I ^ = ej '‘*v 

G G ^ n 


Integrating by parts we obtain 




it 



in G 


* # 

with = &M on O 

<6.19a> can be verified as follows 


lim. p 

^ O J 

ax^ 

1 tffl p p 

= 0 J 

afp^ 

jj 


G* 

3 

i V 

3 

1 M 

3 







Hence = 0 

a^ 

^3 






if if 

We need to obtain relation between ^ and ^ in order to 
the homogenized problem. We define an auxilliary problem 
purpose. For each A c we define as the following. 

8 . t 


(6.18) 


<6. 19) 

<6.19a> 


obta i n 
for this 
F i nd 
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J' = O ¥ V c W 

y * 

► 

and ^ ^ ^ y 

uet = Vu.^ 

Hence from (P-9> we obtain 

' / ~ ^ V V « V 

y 

Pu>^ - h.^ € W 

Let us introduce the functions 
V C}0 = « Cpw.J>Cxy'^^ 

€ - >i - 

2a. = aA‘^2^^ 

Sba. = 

Then from (P-10> we have 

-V ^ 

We can extract subsequences and COg^^} *. t 

w — ► w* in H^CM^ 

€ 

_♦ /< in [L^CM:>]^ 

- € - 

and Vw* = M 


<P-9) 


<P-10> 


< 6 . 20 ) 

< 6.21 > 
< 6 . 22 ) 
<6. 25) 
<6.24> 


An application of Theorem A.II.l CAppendix 11} gives us <6. 25). 

To obtain <6 . 21 ) . <6 . 22) , and <6. 24) we introduce a function V'<iJ) ^ 
« y as follows 

y<^)= Ej)/ = x/«] 

or «V/<x/«5) = ew^Cx/'e? ” ^ 

X + et ^ x / fe ) 


or 
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s i nee v'(^ ) c L 

^(y). v'<x/€) 

«L^(0) 

Therefore w 

► /i. X tn. 

L^CCD 

« 

•* — 

A-- 

2 


-f Vu> »■ h 

— ' — * 

and V w = Ai 

£n L <ra> 



Let. v€ DiG) . Multiplying (6.18) by and (6.20) by we 

obtain 


JJ aj®.7v«.^ + J/ - /J Q3^.7 vP^* 

IW IM IM 


-// v.v<p**) * ;j oi* g ; . //«%; 

IH IM ^ IM 


We note that JJCQ^^ - Qq^^v. VCP^^.^> = O 

IM 

This is as a result of the definitions of the extension 
operators . 

Passing to the limit as e •» O we obtain 


JJ ^ - JJ iny(Qij^>VtJ^* + JJ ?* ^ w*= ©J v*vw 

IM IM 


IM ^3 O 


Integrating by parts and using (6.19) we obtain 

IM 

Since this is true for any v€ D(G) we have 
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5 M = myCQ2^;>.v^ 


Substituting and f* in <6.19) we obtain 


JC«yCS2^:>7^'':>. 7v + J 0^ ~ = 0 f 
G G ^ ^ a 




. a in 


^*-0 

^3 


on 




on n 


The homogenized coefficients are given by the matrix <Q>>j^> . 

Since h is also arbitrary we may choose h = (1,0)^ and (0,1)^ to 
obtain the coefficients in <6.14«) We can now state and prove the 
theorem on the estimates of eigen-values and the final theorem on 
the homogenization problem and conclude the theoritical analysis. 
We shall now prove a lemma which shall be useful in proving the 
homogenization theorem. 

Lomma 6.2 .- Let be an eigen-value of the following problem 


K 


-«? 


ij 9x ,ihc . 
i J 


+ m. ^Cbo^u = u u in O 

y* 


i ,J - 1.2 
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Mhere <?. . are the coefficient defined in <6.14a>. 
Along with the boundary conditions 

on dCi 


9x1 

9tx 


- o 


(P-11 ) 


Let us define w as the solution of the following problem 


-Aw* + ho w* = XI in D* 


9tt> 

- o on 9C1 


<P-12> 


Cwhere w* e 1^(0* >] 

Then the following assertion holds true 

y u irt f/ca:>. (6.25) 

Proof : From. CP-IZ) we obtain 

h'li.o* l«^l l“lc,o 

where Co > o is a constant independent of e. 

From the property of the extension operator P 

we know that |V ?w*|^ ^ < C * 

€ 

With the help of Poincare-Wi rt inger inequality and <6.26> we 


obtain 
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€ 

Hence | ^ ^ 


where > o is independent of €. 

Hence we can extract a subsequence again denoted by € s.t. 

■v f -it i 

P^tT — . w on >r ( 0 > 

It 

We have to show that v satisfies the following equation 




erxjLj 

^ .dx . 

I j 


my 


(&oJ>w 


= 


in O 


9w 


= o 


on ^ Q 


y <P-15> 


Then comparing (P-15> with (P-11) we clearly find that w = 
[from the uniqueness of the solution] 

1 

In such a case we shall have Pu> — 7 w in W <0) and (6.25) is 
proved consequently. 


Let 


£j£ 

3x . 


1.2 


As usual Q shall denote the extension operator which extends the 
function by O in the holes. 

It is clear from (6.26) that ^ C [for some C independent 

of « 3 

* 2 2 

Thus we can extract a subsequence s.t. — 7 in [L (ft) 3 


From <P-12> we have 
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-V.2^* + ^ho^xi^ = in O 


5^ .n = O 


on Xi 


From <P-15) we obtain Vt 7 €D< 0 > 


J ?7 J J ^ 


O 


n 


Cl 


Passing to the limit we obtain 


^ +my<6o) J* w^r) = fjO J* 


Cl 


Cl 


Cl 


•♦ - V. tHyCbo^ w* = BE^ u i n O 

Me now define an auxilliary problem as follows : 

For each A « we define « 14^ as the solution of the 
weak problem 

J ® ^ 

y* 

and Cg^-k.^:> « 1/ Vj; « y 
with a uniqueness condition 


J «A ' ° 


<P-14> 


<6.27) 


<6. 28) 


following 


Let 12/,“? 


Thus from <P-14) we obtain 
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J* ~ ^ 

y 

Let us introduce the follMing functions needed to extend the 
functions defined in the unit cell to whole of £1. 

, = € P Vx € n 

From (P-14) we have 


-V.(Qq^^> = O in £5. 

We can now extract subsequences and CODj^^I s.t. 

7 g* in CO> 

7 inyCQr 7 ^:> in [L^ca>3^ 

and ^g =1$ 


(6.29) 


<6. 30) 
<6. 31) 


(6.32) 

(6.33) 


(6.51) is a direct consequence of Theorem A.II.l [Appendix 113 
To obtain (6.50), (6.52) and (6.35) we first define a function 
y/(Yl as follows: 

¥'<y> * ~ 

or «V/(x/« = ~ ^5" 2£ 

or g <x> = *t. X + 
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since « L^<Y>, v/(x/€> « L^(f5) 

^ z 

Thus ► /i. X in L CC1> 

Consequently we obtain <6.31 > and <6. 33) 

Let x>€ Z)(Q> . Multiplying <6.26) by v g^, <6. 27) by -oPw^ 

subtracting we get [after integrating by parts} 


Jc? ■ ?» 

a 


4c: I ^ ^ e 

Q^T7g^v + X i>o vg^ 

a a 


1-0* ^ € 

I 

O 


“J 2^2^^. vcfSi»*.>v ^ J 

n a 

From the definitions of the operators Q and P we get 


-J 

Q O 


VCPw J>v 


O 


[In the holes ^ a»^cl in the region 


2 ?[^= VCPw^:>3 


Passing to the limit in <6.31) we obtain 


- J«yCQ] 2 ^J>. + JiftyC6o.> V * = JU e 




a 


a 


Q 


£) 


Integrating by parts and using <6.26) we obtain 


and 


6.34) 


and 
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It is clear that the dimension of |^* = I 

for « > O. Here |q* denotes the subspace of whose functions 
are restricted to From the Mini-Max principle of the 
eigen-values we obtain 


Xj = Min 


fi 


Max 
1v€S 


Ci4^u,v^ 1 
I J 


We can express v in terms of Cw^.} as the following 


V = F a .w . where a. is scalar 

» SI s 


L L i, ^ , 

= CA‘^C£ a.w^, E a.CAw.»M>S> ^ CETm, 

J=1 ^ j=i ^ ^ j=t J J J 


v> 


„ f^ox C‘S^v,vJ> 

Hence X < — 

I CB v,vJ> 


II 4. CE^VtvJ> ^ _ 

We now have to prove ~ — — s C. 

V€S, __€ - I 

I CB v,v^ 


<6.55> 


to complete the proof of the theorem. 

We assume the contrary. Then we can extract a subsequence of « 
(again denoted by «> and s. t. 


CE^v^.v^:> = i 


(6.56) 


and -♦ O (6.57) 

It is clear that is bounded in (Q) since is finite 

dimensional. Thus from Rellich's Compactness Theorem 3 a 



133 


subsequence of v* weakly convergent in <fl> and strongly in 
L (Q> . Let us denote the limit by v. 

Thus from <6. 56) we obtain < i and from <6.37> using 

Lemma 6.1 we obtain ( v,xj^= o which is a contradiction. This 
proves <6.55> and hence the lemma. 

We shall now present the final results in the form of two 
homogenize ion theorems. 

Theordm & .2 : 

For any positive integer 121 let X^be the I eigen-value of the 
problem <P-5) and the corresponding eigen-vector. Then 
converges to an eigen-value X^and 3 a sequence of c s.t. 


y in CCD. 


h, ft 

where is a solution of (P-7) corresponding to X 


Proof i Let us denote — in Ci^ i = 1,2 

i 


As usual Q shall denote the extension operator in which the 
function is extended to O in the holes. 

From Theorem 6.1 we conclude that ^ C ^ C where O is a 

constant independent of c] . 

Thus we can extract a sequence of e s.t. 


or®. 7 in L^CCD i = 1,2 

?u® 7 in H^CCD 

X® ► X^ f since from Theorem 1 X® :S 3 
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Now <P-?) can be cast in the following form 


?: + 


I fcO 


r) = \tpT7 Vy> c D<Q) 


(6. >8) 


Integrating by parts we obtain the following operator equation 




<6. 59) 


= o on do* 


Passing to the limit in <6.58) and invoking Lemma 6.1 and Theorem 
A.II.l (Appendix II) we obtain the following 


+ nty<b<?)u^ = X^© J^-u^ in O 


C6.40) 


^ .n =0 


on dO 


We now define an auxilliary problem as follows : 

2 

Determine s.t. for any M € R 


J* Yv = © 


V V € W 


- h. € w Vjjjey 


(P-15) 


Together with a uniqueness condition ° 


Let - V . Then from <P-12) we have 
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-V • Qj 2 ^ = O in O 


<6. 41) 


Let us introduce the following functions needed to 

function defined in the unit cell to the whole of O. 

w = € ?u>. V X € Q 

* « h. - - 

a** = a* 


From <P-12) we have ^ in Q 

Let us introduce a function as follows 


extend 


the 


~ i? V ^ € y 

or = ^ . X + € ^X/€^ = x/€l 

Since € L^iY), yt <x/€> e Z.^ <Q> 

A. 2 

There ^ |5* ^ in L CO) . 

We can now extract subsequences Cw > ca\d <Qti. J s. t. 

4E rtfs 

w 7 w* in cay 

e / 

VW 7 h in [L^cayj^ 

- € / - 


<6.42> 

<6.43> 


-V 

Qj2 




T 


ntyC^^y in [L^cayi 


and tS* = h 


<6.44> 

(6.45) 


Let r»€D<0) . Multiplying (6.59) by vw^and (6.41) by vPu^ and 


substracting we get 
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n o o n 

J = >^tj 2 Jcf*u*^w^ 

o n 


From the definitions of the operators P and Q we get 

J [ 2 ^^^ - 2^3^^. vc?Si»*;>]v = o 

Q 

[Since in the holes 2^^ = = O and in the region 

and Ze^ - VC^i^yi 

Passing to the limit in <6.46> and using <6.42), <6.43>. <6. 

(6.4^> we obtain 

- my • 1}^^ V = O 

Since vs D<.Q) is arbitrary we have 
5*.A = my cZri^:> . 2 

Since At is also arbitrary we may choose M= <0.1)^ and 
obtain the coefficients q. . given in (6.14a}. This proves 


6.46> 


44) . 


to 

the 


theorem. 
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Th»or»m 6.3 : 

Let <M^> be the 

corresponding eigen-vectors 
following hold. 

< i ) X® — 

4b-»^0 


sequences of eigen-values and the 
of the problem <P-7>. Then the 


<i i > 


3 a sequence of again denoted by «) s.t 

Pv^ -* in H^CCD 


Proof : - From Theorem 6.2. <ii> is clear. 

Finally it remains to prove that the limiting point 
X^ of fX^} is the eigen-value of <P-7>. For this it 

suffices to verify that 

<i> there is no eigen-value than for <P-7> <6. 47) 

' cx> Z 

XllJ is an orthogonal basis in L {0>-w.r.t the inner 

product <B^< . > , . ) <6.48) 

For proving (6.47) we argue by contradiction. 

Let us suppose that 3 an eigen-value ij different from X^ with 
eigen-vector which satisfied the homogenized problem <P-7) i.e. 

. 

^ij ^ + my*<t>o)u = p » ri in Q 

i .J = 1.2 
M « V<Q) 
p ^ I 
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and 


1 '^- 


N ^ 1 

)W 1 


We can choose an integer I s.t. ij < k 


1 + 1 


We ^define u> as the solution of the following problem 


-Aw + bo w ~ B XL in O 


9 n 


= O 


on XI 


Then from Lemma 6.2 we have — y m in CCi> Now let us 

€ 

the function v as follows 


V = w 


1 

2 

iv=i 




J = J 


Q 


B \0 XI XL . 

J 


O 


C 

-2 J 


, « € .-C C € 

iB w )B It. u . 

^ J 


o* 


r .€ € € f 


B M XI . - O for j' = 1,2, ... I 


n 


Q 


l+I 


* 2 ^ * 
V « (Q*) 


i4 Cv »v > 


define 


<6. 49) 


Since X 
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i4 Cv »v ^ 2: CB V ,v :> 


lim. r 

c -> O J 


lim 


O 


Q 


w'‘^wr ~ ® J ~ ^ 

Q 


Tharefore Pv^ 7 xt in CCi^ 


.-. A^cv^.v^:> = J V c?u*:>. vc?k>*j> + J bocPo^:>^ 


n 


a 


., lim .€, € € 

Hence ^4 


Cv^.v^^ = J“ V«. Vu + »iy C6oJ> V = /j© J* <^u>u = 


n 


O 


A lim 

Asa.n 


f 3cs^v‘ 


a 


b Pv^ ~ ^ I ci^-u^:>xi = e 
a 


Hence passing to the limit in <6.^0> we obtain 


© 

or fj ^ *^*^^*^^ *® ® contradiction. 

Thus <6.47) is proved. 

To prove (6.48) we proceed in a similar fashion. Define a 
t>^<0> as follows. 


V = w 


J J 


J- = 1,2 ...l 


<6.>0) 


Ai© 


function 


for any 
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Thus we obtain i v,v> 

Since j4^Cv,vJ> < C [ where C > O is some constant 3 


We have CS^v.v^ ► O 

V = O. This proves that C'u^} form a complete orthogonal basis in 
L^<n) w.r.t. the inner product . ) , . ) . 


6.9 NUMERICAL RESULTS : 

We have computed the eigen-frequencies of system 
depicted in Fig. 6.1 for various © Cwhere © is the ratio of the 
dimension of the tube to the dimension of the unit cell!. The 
dimensions of the vessel are 4 x 4 x 4 m and the liquid considered 


is water. 

The finite-element mesh is shown in Fig. 6.9. The 
elements used are the 20 noded 3-D solid elements for the liquid 
and 8-noded isoparametric elements for the free-surface 
discretization. 

The discretization of the representative cell for the 
purpose of computing the homogenization coefficients is given in 
Fig. 6. 4. The variation of the first five frequencies with the 9 
is given in Table 6.1. Five different ratios have been considered 
and for a given 9 the five eigen values are listed in the table. 


6.10 DISCUSSIONS : 

A method has been developed for the computation of the 
eigen-values of the oscillations of an ideal incompressible liquid 




FIG .6.4 2-0 MESH OF THE REPRESENTATIVE 

CELL. 
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FIG • 6.5 


3-0 MESH 
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Table 6.1 

Table showing the first five eigen-values for various ©. 


0.75 0.5 

.1649x10"^ .6561x10 

.1649x10”^ .6361x10 

Py — "5 

.4689x10"^ .1803x10 

.1369x10"^ .5200x10 

.1369x10"^ .5200x10 


Q 


0.25 

0.125 

0.0 

.3724x10"^ 

.3486x10"^ 

.5515x10"'^ 

.5724x10"^ 

.3486x10'^ 

. 3515x10"^ 

.1057x10"^ 

.9904x10"^ 

.9984x10*'^ 

. 5038x10'^ 

.2845x10'^ 

.2868x10"^ 

. 5038x10"^ 

.2845x10"^ 

.2868x10"^ 
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using homogenization techniques. It can be seen that since the 
domain is composed of the translation of the representative cell, 
a single element computation [for stiffness or mass matrices] can 
be used for the global assembly. This brings about a significant 
reduction in the C.P.U. time and the computer memory. 

From the numerical results it is observed that the 
frequencies tend to the frequencies of the system comprising of 
the liquid only and without the tubes as © — ► 0. This is because 
the homogenized coefficients > 1 — ♦ 0 as © -♦ 0- ^ ® 

means that the effect of the tube becomes lesser as the ratio © 
decreases. Also on the other hand when © — + 1 we also 

expect the results to be of much significance since this means 
that the tubes touch each other and the continuity of the region 


is lost. 



CHAPTER VII 


CONCLUSIONS 

7.1 CONCLUSIONS ON THE RESULTS ACHIEVED 

The first problem which has been considered is the 
construction of the capillary surface which in its generality is a 
free-boundary problem. For its solution a shape-optimization 
technique has been developed and applied to some problems. In 
this method the capillary surface problem is reformulated as a 
minimization problem in which the modified energy functional 
[which accounts for the volume constraint] is minimized with 
respect to the domains and a proof for the existence of the 
solution has been given. The major advantage of the method is that 
an accurate and inexpensive technique has been developed to 
compute the gradient of the non-linear energy functional. Hence 
it has become possible to use an efficient optimization algorithm 
based on the gradient for the numerical solution. In the present 
work the ?-noded triangles have been used for the finite-element 
formulation. However it is possible have more accurate 
representation of the surface using higher order isoparametric 
elements . 

The prime task of an adaptive refinement scheme is to 
develop an efficient error indicator. Although the operator 
occurring in the capillary surface problem is strongly monotone, 
the error indicator based on the residual error estimate, does not 
perform satisfactorily due to the fact that the error tends to get 
normalized and regions of high gradients are not detected easily. 
In the work an adaptive refinement scheme based on a new error 
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indicator has been developed which performs satisfactorily. From 
the computational point of view, the efficiency of the method will 
depend on the data structure for the finite element mesh and its 
refining and the method used. 

The problems considered in the work consisted of simple 
geometries [square section! using simple mesh generation methods. 
However, for complex geometries one has to explore other mesh 
generation techniques in order to obtain an efficient one. For 
the non-linear solution the full Newton's Method which involves 
the computation of Jacobian, has been used. Since the Jacobian is 
computed at the element level. the finite-element method is 
ideally suited for the Newton’s Method, although some choice of 
the relaxation parameter can be made from efficiency 
considerations. For the large systems one can use either the 
sky-line solver with block storage or the frontal solvers. 
Although in the work the sky-line solver with block storage scheme 
has been used it may be explored to use frontal solvers in order 
to reduce the C.P.U. time. 

The variational formulation for the problem of the 
linear oscillations of the free-surface becomes difficult due to 
the fact that the differential operator appearing on the boundary 
is of higher order than the differential operator appearing in the 
domain. This problem is re-formulated into a standard eigen-value 
problem with the help of some operators resembling the normal 

2nd mass operators. The difficulty in the computation 
arises in the finite-element modelling of the surface. Inacurrate 
modelling of the surface may give rise to large discretization 
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errors, because they occur throughout the domain and not only near 
the boundary as in the case of planar 2-dimens lonal problems. For 
the eigen-solution several finite-element softwares are available 
which are based on the subspace iteration method or the Lanczos ' 
method. The significant modification which has to be carried out 
in this work is in the assembly of the mass matrix which is fully 
populated unlike the lumped mass matrix. 

The problem of the computation of the eigen-values for 
the liquid oscillations containing large number of rigid bundles 
poses computational difficulties due to the requirement of large 
number of finite-element nodes for its solution. It has been 
solved using the method of homogenization which involves the 
construction of the homogenized operators for the eigen-value 
problem. This is basically a 5-dimensional problem with a 
periodicity in 2-dimensions and requires the construction of the 
homogenized operators both in the liquid domaih as well as at the 
free-surface. The convergence property of the eigen-values and 
the eigen-vectors have been studied in the form of homogenization 

theorems . 


SCOPE OF FURTHER WORK 


There 

aspects of the 
numerical scheme 


is a lot of scope of further work in the 
capillary surface problem. In the 
for the determination of the capillary 


numerical 
work the 
surface 
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was presented for the case in which it is representable as the 
graph of function. The shape-optimization technique so presented 
can be also applied to the generalized 5-dimensional problem. Of 
course in that situation the algorithm will become more complex 
mainly due to the 5-dimensional node movements which requires 
efficient data structures and adequate graphics-support. 

For the dynamic problems, especially the transient 
problems, the number of results is also very small. The theory of 
the dynamic angle of contact is a much debated one and there is 
not agreement on the question whether the dynamic contact angle 
shows a hysteresis phenomenon. 

The finite-element adaptive strategy has been developed 
for the capillary surface problems without singularity. Now the 
method can be extended to study some problems with singularity. 
There is also a good scope of developing a finite-element based 
adaptive strategy for the free-surface flows with surface tension. 
This seems to be a complex task as the surface tension also 
appears in the problem along with the Nav i er-Stokes equation 

The finite-element method which has been used for the 
computations of the eigen-frequencies of oscillations of an ideal 
liquid can be used to study the oscillations of two imiscible 
fluids. It can also be used for studying forced oscillations and 
also the fluid-structure interaction problem taking all phenomena 
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into account viscosity and also behaviour can be studied in 
sufficient details . 



APPEMDIX I 


A. I.l FORMULATION OF THE PROBLEM OF LINEAR OSCILLATIONS OF AN IDEAL 
LIQUID . 

In this section we develop the equations governing the 
free-surface oscillations of an ideal liquid under low gravity 
conditions. The vessel which contains the liquid is assumed to be 
rigid and the liquid to be ideal and incompressible. We also 
assume that the equilibrium configuration of the capillary surface 
to be known a priori. In the sequel we shall follow the notations 
introduced in Chapter V. 

In the liquid region we have the potential flow equation 


= 0 in a 

As regards the boundary condition we have the no-flow condition on 
the wetted part of the vessel surface. 


dn 


O 


on E 


As a kinematic condition we have 


( 2 > 


dip 

9n 


= u on r 


< 5 > 


where xi<x) is the velocity of a point Cwith coordinates x3 on T 
This means that the velocity of the liquid on the free-surface 
matches with the velocity with which the free-surface is displaced 

along the normal. 
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Also since 
cond i t i on 


J 


the volume of the 

“ = ° J Z = 


r r 


liquid 

O 


is conserved we have the 

(4) 


A. I . 2 Dynamic Conditions on the Equilibrium Surface . 


In equation (1-4) the capillary properties (i.e. 
surface-tension (denoted by o) . angle of contact. the curvatures 
of the surface) were not required. These properties are 
manifested in the dynamic conditions on F. We shall apply 
Hamilton’s principle to derive the condition. In the system 
under consideration where the structure is considered to be rigid 
the kinetic energy is provided by the liquid mass and the 
potential energy by the free-surface which behaves like a 
stretched membrane. 

Let p denote the density of the liquid and v the velocity of the 
liquid. Cv = . Let us denote the kinetic energy of the liquid 
for small motion by T. Then 





2 


a 


Let us denote by u the small variation of the potential energy of 
the capillary surface when it is displaced from its equilibrium 
configuration. The expression for xt is given by the following : 



[bon^ 




r 


+ 



Xr? 


2 


dT 
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where 17 is the normal displacement of r [i.e. -a (...) is 

the Beltrami first differential operator [Mysk 1 sC 1 9873 3 . 

Since r is the manifold corresponding to the solution of the 

capillary surface problem. it is known that under certain 

2 4'Oi 

conditions it is sufficiently smooth. In fact it is C I o ^ c*. 

< 13 when the domain of definition of the liquid is convex and 
when its boundary c CMittelmann [19773. Finn [1986b33. Under 

these circumstances the Beltrami Operator is well defined. In our 
work we assume such smoothness on the capillary surface. 


ho = — 


an 

o> "Sri 


_ u2 


on r 


Here n is the gravitational potential and .^2 are the principal 
curvatures at the point. 


X = 


K cos 9 - h 

c 

sin 


on ar 


Here is the angle of contact assumed to be constant along ST. 

and k are the curvatures of the normal cross-sections of the 
surfaces r and E directions (a. 

capillary surface and the vessel surface respectively at the point 
of contact. For further details regarding the derivation of the 
potential energy eypression one may refer to Myskis C19871. 
Considering the variation in T we obtain 


6T = p 








(5) 


n 
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where V = In deriving (?) we have made use of Green’s 

formula along with the boundary condition (2). 

Also for the variation in i/ we have 




J (boT> - 617 + <y J ^ + XTjj Sr) 


< 6 ) 


ar 


While deriving (6) we have used the Green's formula in a manifold 
CMyskisC19873] . 

According to Hamilton's variational principle we have 



Of course we assume that there are not external forces acting on 
the system. The «5^ and Sr) are required to satisfy the isochronism 
condi t i on . i . e . 


Substituting <5) and (6) in (7) we obtain after integrating by 
parts the following equation 
^2 

J [ j[ J C £ * ' ° 

tj r fr 

Including the volume constraint ^Sr) = O with the help of a 
Langragian multiplier y>(,ty we obtain the following modified form 


of <8). 



154 


1 [ ~ - V'<t>]6r? - cx j [ g + Xr>]6r)j = O 




ar 


< 9 ) 


Hence we obtain the following Euler equation from <9) along with 
the boundary condition on dT . 


9<t> 

- a ihan-V^rt) - yf = O in r 


( 10 ) 


^ + Xr? = O on ar (11) 

We divide (10) by a and differentiating <10) and (11) w.r.t time 
variable t we obtain the following equation. 


hou 


Apr, = - 


r 

y in r 

at 


g <»> * X* 


= O on ^ 


( 12 ) 


(15) 


We now assume the harmonic vibrations with frequency and 

substituting in place of if (12) and (12) we obtain the following 

equations : 

ho-u - = K (D^ in r 


<141 


9x1 

9^ 


+ xii = o on ar 


(15) 


[where X = 

along with the volume constraint condition O. Of 

is implied that ^ has to satisfy (1) and the boundary 


course it 
condi t i on 


( 2 ) . 



Appendix II 


In this section we state and prove a theorem which has 
been used quite regularly in chapter VI. 

Theorem A.II.l : Let xitd. (Y ). If we extend it periodically in 0?^ 
have 

E-uC^^ -t- Bty CEtO in Ca> <1) 

-V 

Cwhere E is any extension operator! 

Proof : 

Let = xjtCxye^ V xcCQ?. It is clear that En^ is bounded 

2 

in L (n> as € — > O. Let us choose any t>€D(Q> . Then we have 


f Em^Cx> Cv - V* J> — ► O 

J ^ e. -*■ O 


O 

fwhere v* is a function which assumes in each eY period a constant 
value equal to %> at the centre of the period] 


1 


Evi t> 


my 


Al** 

CE^^;> 



my CExO 

€+0 


1 


V 


I 


Em CxJ> tXrxJ> 


my 


CExO fvCx:> 






we 


Thus (1> is proved. 
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